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1.0  INTRODUCTION 

One  of  the  growing  fields  of  interest  in  the  last  few  years  has 
been  that  of  initiation  of  chemical  chain  reactions  by  the  application 
of  radiation,  with  perhaps  the  greatest  application  in  polymerization 
studies.   To  date,  the  work  that  has  been  done  in  this  field  has  been 
experimental,  with  attempts  to  explain  the  results  in  terms  of  the 
effects  of  overlapping  of  the  clusters  of  radicals  formed  by  the 
radiation  process.   Three  situations  have  been  postulated  as  1)  forma- 
tion of  clusters  of  free  radicals  with  a  separation  distance  large 
compared  to  their  effective  diameter  so  that  most  of  the  radical 
reaction  is  complete  before  their  concentrations  become  uniform;  2) 
formation  of  clusters  of  free  radicals  with   separation  distances 
small  compared  to  their  effective  diameters  so  that  the  initial  distri- 
bution of  radicals  can  be  considered  as  uniform;  3)  formation  of 
clusters  with  a  separation  distance  comparable  in  size  to  their  effec- 
tive diameters  so  that  there  is  some  reaction  between  neighboring 
clusters  in  an  overlapping  region  before  the  complete  randomization 
of  the  radicals.   For  these  cases  the  rates  of  production  of  the  chain 
reaction  product  should  be  a  linear  function  of  the  dose  rate  for  the 
first  situation,  and  range  to  a  function  of  the  square  root  of  the 
dose  rate  for  the  third  situation. 

The  major  purpose  of  this  work  is  to  develop  a  general  computa- 
tional scheme  applicable  to  the  solution  of  the  equation  representing 
the  diffusional  and  kinetic  behavior  of  reacting  species  in  a  single 
cluster,  and  to  apply  this  method  to  the  investigation  of  the  spatial- 


time  behavior  of  a  particular  type  of  chain  reaction.   The  data 
obtained  from  the  history  of  a  single  cluster  cannot  give  actual  num- 
erical knowledge  on  the  effects  of  the  dose  rate  on  overlap  reactions 
but  it  is  possible  to  obtain  estimates  of  the  effects  for  extreme 
cases  of  dose  rates  and  to  indicate  to  some  extent  the  dose  rates  at 
which  overlap  effects  become  important.   Also  of  importance  in  this 
work  is  the  proof  that  overlap  effects  to  some  extent  are  more  impor- 
tant for  chain  reaction  systems  than  for  non-chain  reaction  systems, 
and  the  indication  that  the  use  of  a  modified  form  of  the  prescribed- 
diffusion  hypothesis  should  give  useful  results  even  with  the  complex 
reaction  scheme  considered.   The  results  of  this  work  will  also  form 
the  basis  for  future  computations  where  the  effects  of  spur  overlap 
are  calculated. 


2.0  DEVELOPMENT  OF  THE  THEORY 

2.1   General  Theory  of  Diffusion-Kinetics  in  Radiation  Chemistry 

In  order  to  explain  the  observed  effects  of  radiation  on  media  in 
which  chemical  reactions  are  initiated  by  the  radiation  it  has  been 
found  necessary  to  develop  what  is  commonly  known  as  the  diffusion- 
kinetics  model.   This  model  assumes  that  the  result  of  absorption  of 
radiation  energy  by  a  system  is  the  cause  of  production  of  a  variety 
of  reactive  species  from  the  material  initially  present  with  an  inhomo- 
geneous  spatial  distribution,  these  species  then  undergoing  diffusion 
and  chemical  reaction.   The  over-all  process  is  divided  into  three 
stages.   The  first  stage  is  termed  the  physical  stage  and  consists  of 
the  dissipation  of  the  energy  in  the  system  during  the  time  interval 
of  10    seconds  or  so.   The  second  stage  is  the  physiochemical  stage 
during  which  processes  take  place  that  lead  to  the  establishment  of 
thermal  equilibrium.   This  stage  has  a  duration  of  approximately  10 
seconds  for  aqueous  solutions.   The  third  stage  is  the  chemical  stage 
during  which  the  reactive  species  diffuse  and  react  chemically  until 

chemical  equilibrium  is  reached.   This  stage  usually  has  a  duration 

-8 
of  10   seconds  or  greater. 

The  assumptions  of  the  diffusion-kinetics  model  are  as  follows. 
The  reactive  species  formed  by  the  absorption  of  the  radiation  are 
in  thermal  equilibrium  by  the  time  the  chemical  stage  starts  and  have 
a  specific  spatial  distribution  which  depends  on  the  type  nnd  energy 
of  radiation  used.   These  species  then  diffuse  according  to  macro- 
scopic diffusion  laws  and  react  chemically  according  to  the  same 


rate-laws  that  would  be  obeyed  if  the  species  were  distributed  homogen- 
ously.   Ultimately  chemically  stable  products  are  formed  which  can  be 
determined  by  chemical  analysis. 

In  general  the  deposition  of  the  radiation  energy  in  the  medium 
is  a  highly-localized  process  of  a  statistical  and  complex  nature. 
The  primary  effect  of  this  energy  deposition  is  to  cause  the  loss  or 
gain  of  electrons  by  species  initially  present, or  dissociation,  or  both. 
In  aqueous  solutions  this  leads  to  the  formation  of  the  free  radicals 
H  and  OH  by  dissociation  of  the  water.   The  term  "free  radicals"  is 
restricted  to  molecular  species  in  which  there  is  at  least  one  unpaired 
electron  associated  with  an  atom  of  a  non-metallic  element  whose  valency 
shell  normally  comprises  an  even  number  of  electrons,  all  paired.   The 
presence  of  these  highly-reactive  free  radicals  leads  to  a  system  of 
chemical  reactions  different  from  that  expected  assuming  only  ordinary 
ionization  processes  occur. 

In  order  to  determine  the  initial  spatial  distribution  of  the 
reactive  species,  it  is  necessary  to  know  the  linear  energy  transfer  in 
the  medium  and  to  also  assume  a  value  for  the  amount  of  energy  deposited 
in  each  discrete  interaction  of  the  radiation  with  the  medium.   Once  a 
value  for  the  energy  deposited  is  assumed,  then  the  uncertainty  principle 
limits  the  extent  to  which  we  can  localize  the  wave  packet  associated 
with  the  excitation  produced  by  a  primary  particle.   The  extent  of  the 
uncertainty  is  usually  several  times  larger  than  the  size  of  any  one 
species  so  that  the  excitation  cannot  in  any  reasonable  approximation  be 
considered  localized  in  a  single  molecule.   This  region  is  commonly  termed 
a  "spur".   It  then  becomes  necessary  to  assume  a  form  for  the  spatial 


distribution  of  primary  species  formed  within  this  spur.   The  two  distri- 
butions most  commonly  used  are  Gaussian  and  square,  and  it  has  been  found 
that  either  give  essentially  the  same  results  for  any  diffusion-kinetics 
calculation. 

The  next  consideration  is  the  separation  of  one  spur  from  another. 
Due  to  the  statistical  nature  of  the  interaction  of  the  radiation  with 
the  medium  the  size  of  each  spur  and  the  separation  distance  are  statis- 
tical quantities.   In  order  to  bring  the  complexity  of  the  problem  to  a 
manageable  level  it  is  necessary  to  perform  the  calculations  assuming 
each  spur  is  of  average  size.   It  is  also  necessary  to  perform  the 
calculations  for  two  limiting  cases  only.   Thus  it  is  necessary  to  assume 
that  the  spur  separation  distances  are  large  enough  on  the  average  that 
reaction  between  species  in  different  spurs  is  negligible  compared  to  the 
reactions  within  a  spur,  or  to  assume  that  the  separation  distances  are 
small  enough  that  overlap  between  spurs  occurs  to  such  a  great  extent 
that  a  homogeneous  cylindrical  track  is  formed  before  intraspur  reactions 
are  appreciable.   The  complexity  of  the  computational  problem  is  reduced 
in  either  case  since  only  the  radial  spatial  variable  is  necessary  to 
adequately  describe  the  diffusion  of  the  species  with  time,  as  opposed 
to  the  set  of  three  coordinate  variables  necessary  otherwise.   It  should 
be  noted  that  for  the  case  of  the  cylindrical  tracks,  it  is  also  necessary 
to  assume  that  the  average  distance  between  tracks  is  sufficiently  large 
so  that  track  overlap  effects  are  negligible  and  to  assume  that  the 
tracks  are  sufficiently  long  so  that  track  end  effects  are  negligible. 

One  other  method  that  has  been  used  with  some  success  is  to  assume 


a  form  for  the  intraspur  spatial  concentrations  as  a  function  of  time. 
This  reduces  the  complexity  of  the  problem  such  that  it  is  then  possible 
to  take  into  account  the  statistical  nature  of  the  variations  of  spur 
size  and  separation  distance.   This  method  obtains  its  best  results  when 
the  reaction  rates  are  low  enough  that  negligible  distortion  of  the 
initial  spatial  form  of  the  concentrations  occurs  with  time. 

The  most  serious  criticism  of  the  model  is  that  there  are  so  many 
adjustable  parameters  in  the  solution  that  it  should  be  possible  to  find 
a  set  of  parameters  to  fit  any  result  desired.   This  criticism,  however, 
would  not  be  justified  if  adequate  knowledge  of  the  values  of  the  para- 
meters were  available  from  experiment.   If  all  the  rate  constants  and 
diffusion  coefficients  were  known  the  model  could  be  rigidly  tested  even 
without  accurate  knowledge  of  the  initial  distribution  parameters.   This 
follows  because  for  a  given  medium  and  physical  conditions,  the  distri- 
bution parameters  depend  only  on  the  type  and  energy  of  the  radiation 
used  so  that  the  number  of  experimental  results  to  be  explained  could 
be  much  larger  than  the  number  of  unknown  initial  distribution  parameters. 

In  conclusion  the  test  of  the  validity  of  this  model  depends  on  the 
extent  to  which  it  can  adequately  explain  experimental  results  and  fore- 
cast new  results. 


2.2  Development  of  the  Finite  Difference  Form 
of  the  Diffusion  Kinetics  Equations 

Since  the  development  of  the  partial  differential  equation  describing 
the  time  and  position  behavior  of  reacting  chemical  species  can  be  found 
in  many  references  (1,  2),  it  will  be  sufficient  to  state  in  this  report 
that  the  general  form  of  the  diffusion  kinetics  equation  for  second  order 
reactions  can  be  written  as 

3C±(r,t)      2  N   1   ± 

— ■  V  C  (r,t)  +  I       I     k  C  (r,t)C  (r,t) 

at  j=i  k=i    Jk  J  k 


N 
-C  (r,t)  I      k  C.(r,t),  1=1,2,  ...,  N  (1) 

j=l   1J  J 

where   C.(r_,t)  =  concentration  of  the  i —  species, 

D,  =  diffusion  coefficient  for  the  i —  species, 

i1 

k.,  =  rate  constant  for  the  second  order  reaction  of  species 

j  with  species  k  to  produce  species  i, 
k  .  =  2nd  order  rate  constant  for  the  disappearance  of  species 
i  by  reaction  of  i  with  j , 
V   =  Laplacian  operator, 
N  =  number  of  reacting  species. 
The  term  on  the  left  side  of  Eq.  (1)  is  the  time  rate  of  change  of  the 
concentration  of  species  i.  The  first  term  on  the  right  side  of  Eq.  (1)  rep- 
resents the  change  of  concentration  of  species  i  due  to  diffusion,  the  second 
term  represents  the  creation  of  species  i  due  to  reaction  of  species  j  and  k, 


and  the  third  term  represents  the  disappearance  of 
species  i  by  reaction  of  species  i  and  j . 

The  form  of  the  Laplacian  operator  used  is  that  for  angularly- 
independent  and  axially-independent  cylindrical  tracks,  and  angularly- 
independent  spherical  spurs.   This  form  can  be  written  as 

where      r  =  radial   coordinate, 

a  =  1  for  cylindrical  tracks;   2  for  spherical  spurs. 

For  the  formation  of  free  radicals  by  the  radiolysis  of  liquids  this 
report  will  follow  the  assumption  given  by  Lea  (3)  that  the  initial 
radial  distribution  of  the  free  radicals  is  Gaussian  in  form. 

The  boundary  conditions  are  then  given  as 

Ci(°°,t)  =  0  for  free  radicals,  (3) 

=  C0  for  solutes,  (4) 

30.(0,0 

r =  0  for  radicals  and  solutes,  (5) 

or 

Ci(r,0)   =   <4   for  solutes,  (6) 

No  _r2/r?rRifi 

=  '[2tt(   I)*]*/*   e  for   sPherical  spurs,  (7) 


— ■    °   ,    2      e"r   /t2(Ro)]for   cylindrical  tracks,  (8) 


where     N  =  number  of  free  radicals  in  a  spherical  spur, 


N  /L  =  number  of  free  radicals  per  unit  track  length  in  a 
cylindrical  track, 
C  =  spatially-constant  solute  concentration, 


R  =  empirical  constant, 


Before  choosing  a  finite  difference  scheme  for  the  numerical 
solution  of  Eq.  (1)  it  is  first  necessary  to  modify  the  form  of  Eq.  (1) 
in  order  to  remove  some  of  the  problems  inherent  in  the  numerical 
solution. 

The  first  problem  is  in  the  choice  of  the  time  increments  needed. 
Since,  in  general,  the  concentration  of  any  species  considered  will 
decrease  with  time  due  to  diffusion  and  reaction,  the  reaction  rate 
associated  with  this  species  will  also  decrease.   This  makes  it 
possible  for  the  time  increments  chosen  in  the  finite  difference  scheme 
to  be  increased  with  time  without  decreasing  the  accuracy  obtained  at  each 
time  step.   This  procedure  is  very  desirable  since  it  decreases  the  amount  of 
computational  time  needed  by  a  large  factor. 

There  are  two  methods  of  increasing  the  time  increments.   The  first 
method  is  to  increase  the  time  increment  by  some  arbitrary  factor  several 
times  throughout  the  computation.   The  points  at  which  the  time  increments 
are  increased  also  must  be  chosen,  based  upon  some  suitable  criterion. 
It  is  difficult,  however,  to  choose  a  criterion  that  satisfies  all 
reaction  possibilities,  so  that  it  is  necessary  to  limit  the  criterion 
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to  the  particular  reaction  system  chosen.   Even  with  this  limitation,  a 
suitable  criterion  can  only  be  found  by  a  trial  and  error  process. 

The  second  method  for  increasing  the  time  increments  is  analytical 
in  nature.   One  attempts  to  find  some  suitable  relationship  between  time 
and  another  variable  such  that  constant  increments  in  the  second 
variable  correspond  to  increasing  increments  in  time.   This  relation  is 
substituted  into  the  original  partial  differential  equation  and  a  finite 
difference  scheme  used  with  the  second  variable.   A  relationship  which 
is  suitable  for  the  reaction  system  considered  in  this  report  is  that 
given  by  Dyne  (4)  as 

t  =  t0e\   o<.x<°°,   t^t^eo,  (9) 

where     t  =  time, 

X  =  dimensionless,  independent  variable, 
t  =  scaling  factor  corresponding  to  time. 
The  problem  now  becomes  one  of  choosing  the  proper  t  and  X -increment 
such  that  the  desired  accuracy  is  obtained  for  the  solution.   This  is 
once  again  a  trial  and  error  procedure. 

The  second  problem  that  must  be  resolved  is  that  of  keeping  pace 
with  the  expansion  of  the  free  radicals  in  the  radial  direction  as  time 
increases.   An  analytical  solution  to  this  problem  is  given  by  Dyne  (4) 
as  the  introduction  of  a  new  variable,  p,  replacing  the  radial  variable 
in  such  a  fashion  that  a  constant  range  in  p  corresponds  to  a  time- 
increasing  range  in  r,  the  radial  variable.   The  range  in  r  must  increase 
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at  the  same  rate  that  the  radicals  diffuse  outward  in  order  that  signif- 
icant information  is  not  lost  from  the  problem.   It  can  be  shown  (see 
Appendix  A.l)  that  the  proper  form  for  p  for  Gaussian  spurs  or  tracks  is 
given  by 

p  =  r/vT7F  ,    t0<t<oo  .  (10) 

The  constant,  t0,  is  now  given  as  (see  Appendix  A.l) 

tQ  =  R0/(2D).  (11) 

The  substitution  of  Eq.  (9)  and  Eq.  (10)  into  Eq.  (1)  will  give  the 
required  form  for  the  diffusion  kinetics  equations  provided  that  all  the 
diffusion  constants  and  initial  Gaussian  widths  at  half-maximum  are  the 
same.   If  this  is  not  the  case,  then  the  problem  arises  of  how  to  best 
follow  all  of  the  diffusing  radical  species  with  time.   One  method  is  to 
define  a  set  of  p  ,  one  for  each  diffusing  free  radical  species,  and  to 
use  an  interpolation  formula  in  computing  the  reaction  terms.   In  this 
manner,  the  errors  introduced  in  computation  of  the  diffusion  portion  of 
the  kinetics  equation  for  each  species  will  be  minimized  as  much  as 
possible,  since  each  species  will  be  followed  independently  in  position. 
The  error  introduced  by  the  use  of  an  interpolation  scheme  in  order  to 
compute  the  reaction  terms  will  also  be  minimal  due  to  the  well-behaved 
spatial  forms  for  the  concentrations.   This  method  will  tend  to  increase 
the  problem  computation  time  due  to  the  need  for  interpolation,  but  at 
the  same  time  the  accuracy  will  be  significantly  increased. 
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If  we  now  redefine  Eqs.    (9),    (10),    and    (11) 


t  E    t^te*-!]     ;        o<X<^    ,  (12) 


p<    E  r:  (13) 

1    /y^+t)   • 

i  =    (Ro)2  CIA) 

Substituting  Eq.  (13)  into  Eq.  (1),  we  have  for  the  various  terms 


aci(r,t)   aci(Pi,t)  3Pi     1    aci(pi,t) 

3r        9p     3r   "  /D^t^t)    3Pi 


(15) 


32Ci(r,t)  1  3      3Pi   3Ci(pi,t)  1  32Ci(pi,t)  (l6) 

3r2  /Di(ti+t)      3Pi    3r  3p±  ^(tj+t)        3p±2 


3Ci(r,t)         3Ci(Pi>t)         3Ci(Pi,t)    3p± 
3t  3t  3pi  3t 

3Ci(p1,t)  r 3Ci(Pi>t) 

3t  2(ti+t)/Di(ti+t)  '       3Pi 

9Cj(pj,t)  Pi  9Cj(p1,t) 

3t  2(tj+t)  3Pi  <17> 
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Eq.  (1)  becomes 


3C.(p.,t)    92C. 
(ti+t) : =  r  H 


3t 


9p 


pi 


+ 


9C, 


■i-3pi 


+  (tj+t) 


1  hi 

;.  =  1  m=l 


C£(Pi,t)Cm(Pi,t)  - 


CiCPi.t)  i  kuCa(Pi,t) 
£  =  1 


(18) 


Substituting  Eq.  (12)  into  Eq.  (18)  gives 


ax 


1+e 


xr?i 


ti 

-o 


-  1 


32Ci 
— 5~  + 
3P. 


Pi   a 


3pi 


1  X 


toe 


N    I  N  _ 

I   £  k£mC^Pi»x)Cm(pi,x)  "  CiCpi.x)  [  kl£C£(p  ,X) 
£=1  m=l  1=1  1 


(19) 


Eq.  (19)  is  in  the  proper  form  such  that  a  constant  increment  in  x  will 
give  an  exponentially- increasing  increment  in  time,  and  a  constant 
increment  in  p  will  give  time-increasing  position  increments  such  that 
the  same  portion  of  the  Gaussian  concentration  for  the  free  radicals 
will  be  utilized  at  each  time  step. 

In  order  to  solve  the  set  of  nonlinear  partial  differential  equations 
given  by  Eq.  (19),  it  is  Convenient  to  use  some  form  of  finite  dlffi 
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for  the  various  terms.   The  choice  of  the  difference  scheme  used  is 
determined  by  experience  as  there  are  many  difference  schemes  available. 
The  scheme  used  in  this  report  is  that  of  central  differences  in  P  and 
forward  differences  in  X. 

If  we  let  Ap  be  the  finite  increment  chosen  along  the  P  axis,  and  j 
be  the  index  associated  with  this  axis,  and  if  we  let  Ax  be  the  finite 
increment  chosen  along  the  x  axis,  and  k  be  the  index  associated  with 
this  axis,  the  derivatives  in  Eq.  (19)  become 

9Cj(p  j,X)  m   Cj(j,k+1)  -  Cj(j,k)  f  (20) 

3X  AX 

3Cj(Pi,x)  =  Cj(j+l,k)  -  Cj(j-l,k)  ^  (21) 

dPi  2Ap 

92Cj(Pi,x)  =  Cj(j+l,k)  -  2Cj(j,k)  +  Cj(j-l,k)   .  (22) 

3Pi2  (Ap)2 

where     Cj(j,k)  =  concentration  of  species  i  at  the  location  p  =  jAp 

and  X  =  kAX. 

It  is  to  be  noted  that  the  truncation  error  caused  by  the  approx- 
imations to  the  derivatives  with  respect  to  p^  is  of  order  (Ap)2,  and 
the  truncation  error  caused  by  the  approximation  to  the  derivative  with 
respect  to  X  is  of  order  (AX).   There  are  several  reasons  why  a  better 
approximation  scheme  is  not  chosen  for  the  derivative  with  respect  to  X, 
the  first  being  that  any  other  scheme  makes  the  starting  of  the  problem 
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solution  difficult  due  to  the  lack  of  sufficient  initial  data.   For 
instance,  if  a  form  for  Eq.  (20)  were  chosen  similar  to  that  for  Eq.  (21), 
it  would  still  be  necessary  for  the  first  step  in  X  to  be  computed  using 
the  form  given  by  Eq.  (20)  in  order  to  start  the  solution.   The  second 
reason  for  using  the  given  form  is  that  of  limited  computer  core  storage. 
Any  other  form  for  the  derivative  will  require  that  additional  values  of 
C.  be  stored  in  the  computer  as  the  problem  is  solved.   The  third  reason 
for  the  form  chosen  is  concerned  with  the  stability  of  the  solution  and 
is  discussed  in  the  section  of  this  report  on  stability. 

Before  substitution  of  Eqs.  (20)  through  (22)  is  performed  there  is 
a  problem  that  exists  with  Eq.  (19)  and  that  is  the  singularity  at  Pj_  =  0. 
Using  L'Hospital's  rule  we  have 


Lim 
P-i+o 


1  3Ci(Pi,x) 


3Pi 


Lim 
P-f^o 


H2 


3zci(Pi,xy 


3pi' 


32Ci(Pi,x) 
3Pi< 


iPi=o 


(23) 


This  limit  is  allowable  since 


Lim   Ci(pj,x)   _ 


=  0  due  to  the 


p-i+o     3pi 

symmetry  of  the  concentrations  around  r  =  0  (hence  around  p^  =  0,  since 
P^  is  a  linear  function  of  r) . 

Eqs.  (20)  through  (23)  allow  Eq.  (19)  to  be  written  in  finite 
difference  form  as 


For  p-^o: 
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Ci(j,k+1)  = 


Ci(j+l,k) 


AX 


+  iM.  +    aAx 


(Ap)2    A    2j(Ap)2 


+  Ci(j-l,k) 


AX     JAx    ctAx 
(ApT    4    2j(Ap)2 


2AX 


(Ap)2   X 


C,  (j,k) 


1+e 


-kAX 


+  CjCj.k)  +  AXt^ekAx 


N    £ 
A=l  m=l 


-  C±(j,k)  I   ki^qCPi.k) 
£=1 


+  C.(j,k) 


(24) 


ForPj^  =  o: 


CjCo.k+l)  = 


2(l+a)AX 


1+e 


-kAX 


■Eft-  1 

to 


C  (l,k)  -  C  (o,k) 
[Ap]2 


+  tiekAX 


AX 


~N    £   .  N  _ 

I        I   k£mC$>»k)Cm(°»k>  "  Ci(0»k)  I  ki£C£(°>k) 
£=1  m=l  A=l 


+  (^(o.k)  . 


(25) 


Inspection  of  Eq.  (24)  shows  that  another  item  must  be  considered, 
and  that  is  the  reaction  terms  of  the  form  C  (p-,k)C  (p-,k).   In  general, 


m 
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unless  the  diffusion  constants  and  initial  Gaussian  widths  at  half- 
maximum  are  the  same  for  species  %   and  m  as  those  for  species  i,  no 

values  for  C„(p.,k)  and  C  (p,,k)  will  be  obtainable  at  p.  due  to  the 
I      1'        mi  1 

nature  of  the  finite  difference  scheme.   For  this  reason  it  is  necessary 
to  introduce  some  sort  of  interpolation  scheme  in  order  to  compute  the 
reaction  terms. 

Writing  Eq.  (13)  for  species  i  and  I   we  have 


Pi  = 


A 


D±(tjft) 


(26) 


PZ   = 


A 


Df(t04t) 


(27) 


Combining  Eqs.  (26)  and  (27)  to  eliminate  r 


P£  =  Pi 


D,(tVt) 


-l1^ 


(28) 


This  last  equation  then  essentially  represents  a  spatial  coordinate 

transformation  from  the  i   system  to  the  I        system.   Note  that  in 

general  —  will  not  be  an  integer  so  that  interpolation  is  necessarv. 
Ap 
Since  the  spatial  distributions  encountered  in  this  report  are 

smoothly-varying  functions,  a  simple  second-order  interpolation  scheme 

is  used.   For  those  cases  where  the  interpolated  quantity  falls  outside 

the  range  considered  in  the  numerical  solution,  the  value  is  tnken  to 

be  zero. 
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The  final  detail  considered  in  this  section  of  the  report  is  that 
of  how  best  to  follow  the  solute  concentrations  in  time  and  space.   The 
requirements  here  are  different  from  that  for  the  radicals  due  to  two 
facts:   (1)  solute  concentration  changes  are  primarily  due  to  reaction 
rather  than  diffusion,  and  (2)  the  solute  concentration  does  not  go  to 
zero  far  from  the  center  of  the  spur  or  track.   In  order  to  insure  that 
the  solute  concentrations  are  well  followed  in  time  and  space,  it  is 
necessary  that  a  spatially-increasing  coordinate  system  be  chosen  for 
each  solute.   The  system  used  in  this  report  is 


PH  m   j ; (29) 

where  the  i   superscript  here  refers  to  species  i.   Species  i  is  that 
species  reacting  with  solute  %   which  has  the  largest  diffusion 
coefficient. 

If  we  use  Eq.  (29)  and  perform  a  derivation  similar  to  the  one 
done  previously,  we  obtain  the  following  finite  difference  equations  for 
the  solute  concentrations. 
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For      £>o: 


C£(j,k+1)   =  CA(j,k)  + 


Dt  ax 


D  -kAX 

1    1+e 


rii.r 

1.1 


C„(j+l,k) 


(Ap) 


J 

j 

-  + 
4 

— 
a 

2j(Ap)2 

+  C£(j-l,k) 


1      j 


(Ap)2        A        2j(Ap)2 


(Ap)2 


-c.a.k) 


+AXt^e 


kAX 


N      m     . 
m=l  n=l 


m=l 


(30) 


For      I  =  o: 


C£(o,k+l) 


C£(o,k)  + 


D£  2(l+a)AX 


'i  -kAX 

1+e 


4 


c^a.k)  -  c£(o,k) 

(Ap)2 


+  AXtie 


kAX 


N        m 

V        V  k1  C    (o,k)C    (o,k) 
m=i   n=l 


N 
Ct(o,k)    I  k£mCm(o,k) 
m=l 


(31) 
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Note  that  the  only  difference  between  Eqs .  (30)  and  (31)  and  Eqs . 

(24)  and  (25)  is  the  addition  of  a  factor  — —  .   In  practice  it  is 

i 

usually  sufficient  to  set  D.  =  D  so  that  Eqs.  (30)  and  (31)  are 

identical  with  Eqs.  (24)  and  (25)  previously  derived.   This  simplifies 

the  programming  somewhat,  and  will  give  sufficient  accuracy  unless  D. 
and  D  are  quite  different  in  value 
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2.3  Development  of  the  Finite  Difference  Form 
of  the  Reaction  Integrals 

The  reaction  integrals  give  the  total  amount  of  a  particular  species 
formed  or  destroyed  by  a  particular  reaction  and  are  important  in 
calculating  yields  and  checking  the  stability  of  the  numerical  solution. 
These  integrals  are  given  as 


t 

r 


V°  =  j 


dt 


dV  k..C.(r,t)C.(r,t)  ,  (32) 


and 


N  (t)  =  JdV  C.(r,t)  ,  (33) 

0 

in  which  L   (t)  »   the  number  of  the  i —  radicals  created  bv 

reaction  of  species  i  with  species  j 

in  a  single  spur  up  to  time  t.   L   is  the  number  of 

i —  radicals  created  with  k , .  replaced  by  k.,  , 


N  (t)  ■  the  number  of  radicals  of  species  i  present  at  time  t 
in  a  single  spur,  or,  the  number  of  radicals  per  unit 
track  length  present  at  time  t  in  a  single  track; 
dV  =  differential  volume, 
■  4:tr2dr  for  spurs, 
c  2irrdr  for  tracks. 


22 


As  a  check  on  the  solution  obtained  a  material  balance  gives  the 
following  relationship  which  must  hold  at  every  time  t  : 


"jW^W-hijW  +  ^j^t). 


(34) 


Deviations  from  this  equality  are  inherent  in  any  numerical  solution  and 
the  extent  of  deviation  can  only  be  used  in  a  qualitative  manner  based 
upon  experience  to  form  an  opinion  on  the  actual  accuracy  of  solution. 

In  order  to  put  Eqs.  (32)  and  (33)  into  a  form  consistent  with  the 
calculational  procedure  for  the  solution  of  the  diffusion-kinetics 
equations  Eqs.  (12)  and  (13)  are  used  in  the  form 


dt  = 


o 


C~e  dX  ' 


(35) 


dr  = 


_ dp.   . 

M).  ft1  +  t^eX-ir 


(36) 


Substituting  Eqs.  (35)  and  (36)  into  Eqs.  (32)  and  (33)  gives 


For  spurs : 


Vx> 


^(D.)3^.  }dxe: 


tj  +  tl(eX-l)J"! 


•  [dPiPiCi<Pi>x>Ci<Pi'x> 

4 


(37) 
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N^x)  =  Att  D1 


t1  +  tl(eX-l) 

0     o 


3/2 


0 


dp.p?C  (p.,x)   , 
1  1  i  1 


(38) 


For  tracks: 


Li:j(x)  -  2*0^. 


dXeX 


ti  +  t1(eX-l) 

00 


dPiPici(Pi,x)cj(Pi,x)  , 


(39) 


BO 

Nt(x)  =  2KD.  ti  +  tl(eX-l)j  Jdp.p.C^p., 


x) 


(40) 


In  order  to  numerically  evaluate  these  integrals  they  are  written 
in  the  following  form 


Lij 


M  7 

■      I  AX  JdpiFiCi(pi,mAX)Cj(pi,mAX)kijtlemAX      ,  (41) 


m=l 


»i<X>    "      jdpiFiC1(p.,x) 


(42) 


where  M  =  x/AX,    an   integer, 

Fi   =   Ati   D1]tJ  +   t^(emAX-l)"|  for   spurs, 


3/2 


«i[<: 


2wD4lt  +  t1(etnAX-l)~|   for  tracks   . 


) 1   for  tra< 
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In  order  to  evaluate  the  volume  integrals,  a  second-order  numerical 
scheme  is  used  (see  explanation  of  "FATES"  subroutine  elsewhere  in  this 
report)  which  gives  a  weighting  factor  WNR(£)  such  that  Eqs.  (41)  and 
(42)  can  be  written  as 

-  M  A 

L..(X)  =  I   AX  I   FiCi(i,m)C .(P  ,m)k  ±  t*e  *WNRU)   ,  (43) 

J      m=l   *=1         J       J 


N.(x)  =   7  F,C.(£,M)WNRU)   .  (44) 

1      lk   X   X 


Once  again  we  notice  that  in  general  the  term  C.(p  ,m)  will  not  have 
tabulated  values  so  that  interpolation  is  necessary.   In  this  case, 
however,  the  quantity  C. (£,m)C. (p . ,m)  has  previously  been  calculated  in 
the  solution  of  the  diffusion-kinetics  equations  so  that  it  is  only 
necessary  to  use  this  value. 
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2.4  Stability  of  the  Finite  Difference  Form 
of  the  Diffusion-Kinetics  Equations 

Since  the  exact  analysis  of  the  stability  of  the  given  set  of 
nonlinear  partial  differential  equations  used  in  this  report  is  extremely 
difficult  if  not  impossible,  stability  criteria  will  be  determined  only 
for  the  case  where  there  are  no  reaction  terms.   The  possible  effects  of 
the  reaction  terms  will  then  be  discussed  in  a  qualitive  fashion. 

Writing  Eqs.  (24)  and  (25)  in  matrix  form  without  reaction  terms, 


where 


C.(k+1)  =  B^.Ck) 


(45) 


C.(k)  = 


C.(o,k) 

(^(1,10 

c^U.k) 


(46) 
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-2(l+a)Ax 


(Ap) 


1+e 


-kAX 


AX 


(AP) 


n   +  1 


V 


1+e 


-kAX 


rt 


ti 


-i 


2(l+a)AX 


(AP)2 


1+e 


-kAX 


£-1 


-2AX/(Ap) 


AX 


1+e 


-kAX 


r>i 


+  1 


L(AP 


_(AP) 


+  1  + 


2(AP) 


^-1 


1+e 


-kAX 


f-1 


^LCApT7""  2  "  A(Ap)tJ 


-2AX/(Ap) 


1+e 


-kAX 


i_o  -J 


1+e 


-kAX 


r&-r 


<«) 


Now  a  necessary  and  sufficient  condition  for  stability  of  the  given 
set  of  equations  is  that  all  of  the  eigenvalues,  x>  of  the  matrix  B^ 
satisfy  the  relationship  |x|<l  •   In  order  to  obtain  bounds  on  the 
eigenvalues  we  can  use  the  Gerschgorin  theorem  which  states  that  the 
largest  eigenvalue  is  equal  to  or  less  than  the  maximum  value  of  the  sum 


of  the  magnitudes  of  the  elements  in  any  row.   That  is,  if  B  =  (b.,)., 
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then 


[X   ]  .<  max  [  [b1k].   . 
max  1  -  .   .    J K  1 


(48) 


Applying  this  theorem  to  the  first  row  we  have 


AX 


(49) 


(Ap)2   2(l+a) 


1+e 


-kAX 


Lro   _j 
th 


-1 


Applying  this  theorem  to  the  j   row  we  have 


AX 


(50) 


(Ap)2  "  2 


1+e 


-kAX 


4 


-i 


-i 


In  order  to  obtain  Eqs.  (49)  and  (50)  it  was  necessary  to  assume  that 


ti  >  tj  . 
o  —  o 


(51) 


It  is  also  easily  seen  that  the  most  limiting  requirement  is  that  given 
by  Eq.  (49)  so  that  this  is  the  criterion  used. 

Now  it  is  obvious  that  the  smallest  value  for  the  terra  on  the  right- 
hand-side  of  Eq.  (49)  is  obtained  when  the  value  of  k  is  infinite.  This 
gives  the  criterion  that  must  be  met  as 
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AX    <    I  . 

(Ap)2"2(l+0)  (52) 


For  the  case  where  a  =  2,  this  requirement  may  be  too  restrictive.   This 

can  be  seen  by  noting  that  in  this  case,  the  magnitude  of  the  first  term 

in  the  second  row  of  matrix  B  will  become  =  A^M  which  will  be  of  the  order 

0.01  or  so.   If  we  take  this  value  to  be  identically  zero  in  comparison 

to  the  magnitudes  of  the  rest  of  the  quantities  in  the  matrix,  the  stability 

criterion  becomes 

_^X_<I  m  (53) 

(Ap)2  ~3 

The  criterion  given  by  Eq.  (53)  is  probably  more  realistic  than  that 
given  by  Eq.  (52)  for  spur  calculations. 

It  is  also  of  interest  to  compare  the  criteria  obtained  here  with 
that  obtained  by  Kupperraann  (5)  for  the  case  where  finite  differences 
were  used  directly  with  the  original  diffusion-kinetics  equations  in  time, 
t,  and  space,  r.   From  Eqs.  (12)  and  (13)  we  obtain  the  following 
transformation 


Didt  m  tl*X  dX    .  (54) 

i  _  t-1  +  t-UX  (An.s2 


(dr)2   tj  -  ti  +  t^eA  (dPi) 
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If  we  let  i  =  1  and  replace  the  derivatives  with  finite  difference 
approximations  we  have 

DlAt   „   AX  (55) 

(Ar)2    (Ap)2 

Since  the  equality  given  by  Eq.  (55)  is  approximately  true,  it  follows 

that  the  criterion  given  by  Eqs.  (52)  and  (53)  should  be  very  nearly  the 

same  as  that  given  by  Kuppermann  for  the  quantity  D^At   .   This  would 

(Ar)2 
probably  not  be  the  case,  however,  if  the  type  of  finite  difference 

scheme  used  by  Kuppermann  was  significantly  different  than  that  used  here. 

It  can  be  seen  by  referring  to  Kuppermann1 s  paper  that  the  set  of 
criteria  are  in  fact  identical  for  the  case  where  i  =  1,  thus  further 
substantiating  the  validity  of  using  Eq.  (53)  for  spur  calculations. 

Before  proceeding    to  a  general  discussion  of  stability  with 
reaction  terms  present,  it  is  necessary  to  show  the  change  in  the 
stability  criteria  obtained  if  the  following  more  accurate  approximation 
to  the  time  derivative  were  used  : 

3Ci(p1,x)  a   C.(j,k+1)  -  Ct(j,k  1)   >  (56) 

3X  2AX 

This  approximation  will  cause  two  changes  in  Eqs.  (24)  and  (25).  All  of 
the  terms  on  the  right  side  of  these  equations  will  be  multiplied  by  two 
except  for  the  Ci(j,k)  term,  which  will  be  replaced  by  Ci(j,k-1).   If  we 
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now  assume  that  the  following  relation  holds  everywhere 


Ci(j,k-i)  =  ejkCi(j,k)    , 


where  B   is  a  number  very  close  to  one,  we  have  a  set  of  equations  of 
jk 

the  form  given  by  Eq.  (45)  and  can  obtain  the  stability  criterion  as 


(57) 


AX  <   1  (58) 

(Ap)2  ~4<l+o) 


Comparison  of  Eq.  (58)  with  Eq.  (52)  shows  that  the  more  accurate 
approximation  effectively  doubles  the  number  of  time  mesh  points  needed 
for  a  given  spatial  increment  in  order  for  the  stability  criterion  to  be 
satisfied.   This  is  a  major  reason  for  using  the  less  accurate  finite 
difference  approximation. 

Concerning  the  influence  of  the  reaction  terms  on  the  stability  of 
the  solution,  a  few  generalities  can  be  obtained  if  the  reaction  terms  are 
written  in  the  form 


C1(J,W61(J.k)  =  AXtiekA*  I        I   ^(p  k)C  (p  k) 

£=1  m=l 


N 
-  C.(j,k)  I   ki£C£(pifk)    ,  (59) 

£=1 


where  it  is  assumed  that  it  possible  somehow  to  determine  G-(j,k) 
The  stability  criterion  obtained  for  this  case  is  given  as: 
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Ay    2-max|G.  (j  ,k) 


(Ap)2 


4(1+ a) 


1+e 


-kAX 


tl 
o 


(60) 


From  Eq.  (60)  it  is  easy  to  see  that  in  order  to  keep  the  original 
stability  criterion  valid,  max|G^(j,k)|  must  remain  small  in  comparison 
to  unity.   Physically  this  means  that  the  change  in  the  spatial  concen- 
tration profile  from  one  time  step  to  the  next  caused  by  reaction  must 
be  small,  at  least  in  the  regions  where  the  concentrations  of  the  various 
species  are  large  so  that  the  possible  growth  of  errors  can  cause  a 
significant  error  in  the  final  answer.   The  maximum  allowable  value  for 

G  (j,k)  is  determined  once  a  value  is  chosen  for  — ^ —  .   In  actual 

1  "  (Ap)2 

practice,  however,  it  is  usually  sufficient  to  set  a  value  for  Ap  and 

determine  Ax  from  Eq.  (60)  for  the  case  of  j  =  k  =  0  (that  is,  from  the 

given  initial  conditions  in  the  center  of  the  spur  or  track  where  the 

reaction  rate  is  greatest).   This  procedure  will  be  sufficient  generally 

for  the  entire  solution  for  the  case  where  the  reaction  rates  continually 

cause  decreases  in  the  values  of  the  concentrations  of  all  the  species 

since  the  absolute  value  of  G^(j,k)  will  generally  decrease.   Unfortunately, 

the  opposite  is  also  generally  true  for  the  case  where  the  reaction  rates 

cause  the  concentration  of  one  or  more  species  to  increase.   For  this 

case,  it  may  become  manditory  to  change  the  value  of  A\  several  times 
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during  the  course  of  the  solution  in  order  to  keep  the  errors  within  a 
desired  level.   The  decision  about  what  type  of  computation  scheme  to  use 
for  any  particular  reaction  scheme  can  only  be  reached  after  much  consid- 
eration and  experimentation. 

For  the  reaction  scheme  considered  in  this  report,  it  is  sufficient 
to  determine  a  value  for  Ax  from  Eq.  (60)  for  j  =  k  =  0  and  to  leave  this 
value  unchanged  throughout  the  entire  calculation. 
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2.5  The  Reaction  Scheme  Considered  and 
Development  of  the  Equations  Describing  the 
Rate  of  Product  Formation  as  a  Function  of 
Dose  Rate 

The  reaction  mechanism  considered  in  this  report  is  given  by  the 
following  set  of  equations 

R2  +  Radiation >  R  +  R  ,  (61) 

kll 
R  +  R — >  R   ,  (62) 

2 

km 

r  +  S >  RS  ,  (63) 

R2  +  RS      >  P  +  R  (64) 

where     R  =  a  free  radical, 
R2  ■  a  solute, 
S  =  a  solute, 
P  =  a  final  product, 
RS  =  an  intermediate  product, 
kj.  =  rate  constant  for  reaction  (62), 
k.   =  rate  constant  for  reaction  (63), 
k   ■  rate  constant  for  reaction  (64). 
The  first  equation  represents  the  formation  of  free  radicals  by  radiolysis 
of  the  solute,  R   to  give  the  initial  Gaussian  distribution  of  free 
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radicals  used  in  the  calculations.   The  consequences  of  this  reaction 
appear  in  the  initial  conditions.   The  second  equation  describes  the 
radical  recombination,  or  chain-breaking,  reaction.   The  last  two 
equations  describe  the  chain  propagation  step. 

Letting  a  bracket  denote  the  concentration  of  each  of  the  species, 
the  diffusion  kinetic  equations  to  be  solved  are 


-^!  =  DRV2[R]    -  kn[R]2   -  kllf[R][S]   +  k23[R2][RS]  ,  (65) 

3 1 


3[R  ] 

2-  =  DR  V2[R2]    +  %kn[R]2    -  k23[R2][RS]    , 


=  D  V2[S]    -  k..  [R][S]    , 

at     s        1H 


(66) 


(67) 


3[RS] 


=  D  V2[RS]  +  k1J+[R][S]  -  k23[R2][RSl  ,  (68) 


— -=  DPV2[P]  +  k,,[Rj[RS]  .  (69) 

dt     p         23  2 


In  order  to  describe  the  forraation-of-product  rate  as  a  function  of 

dose  rate  for  the  reaction  scheme  considered  it  is  necessary  to  integrate 

over  spatial  dependence  of  Eqs.  (65)  through  (69).  Applying  the  spherical 
volume  integral  to  the  diffusion  term  in  Eq.  (65)  we  have 
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dVD  V2[R]   =   4ttD 


fdr  r«!!!U 


8ttD 


3r< 


dr   r. 


3[R] 


3r 


(70) 


Integrating   the   first   integral   on  the  right   side  of   Eq.    (70)    gives 


'dr  r2ii!ii  =     r2ii*l 


Zrz  3r 


OO  00 


-2 


frdri[!li 
3r 


(71) 


Substituting  Eq.  (71)  into  (70)  gives 


00 

fdVD  V2  [R]  =  4ttD     r2 


9[R] 

3r 


(72) 


For  the  cases  of  either  Gaussian  conditions  or  solutes  it  is  seen 


that 


Lira        „3  [R]  Lira        ,3  fe  ] 

r^ — ~ —  ■      r^ 

r-*o     3r     r-*30     3r 


=  0 


(73) 


Therefore  Eq.  (72)  becomes 


dVD  V2[p  ]  =  0. 

R 
V 


(74) 


For  each  of  the  reaction  terms  we  make  the  following  definition 
for  brace  notation 


{[R](S]}  = 


dV(R](s]  . 


(75) 
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This  gives  for  the  set  of  differential  equations 


d{[R]} 


dt 


=  -  kn{[R]2}  -  kllt{[R][S]}  +  k23{fR2][RSP    ,  (76) 


^-^2—=  hkn{[R}2)   -  k      {[R2][RS]}    ,  (77) 


2  3 
dt 


d{[S]} 

1   =  -  kllt{[R][S]}    ,  (78) 


dt 


d{  [RS]} 

=  k1If{[R][S]}   -  k23{[R2][RS]}    ,  (79) 


dt 


^^  =  k23{[R2][RS]}  .  (80) 

dt 


For  the  case  of  low  dose  rates  we  will  assume  that  the  chain  length 
is  large  so  that  a  near  steady  state  position  will  be  reached  with  respect 
to  the  intermediate  product,  RS.   This  means  that  the  radical  recombination 
is  negligible  and  a  balance  is  reached  where  the  total  number  of  radicals 
and  intermediate  product,  RS,  remain  nearly  unchanged  with  time.   For 
this  case  we  set 


d{[RS]} 

-  0  =  k      {  [R][S]}   -  k      {  [R    1[RS]}    .  (81) 

dt  lk  23        2 
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Substituting  into  Eq.  (69)  gives 


k  {[R][S]}  .  (82) 

dt      1£* 


If  we  now  assume  that  the  total  amount  of  solute,  S,  remains 
approximately  spatially  constant  at  a  concentration  [S]   we  have 


^-L-k^tsyiR)}  .  (83) 


This  gives  the  rate  of  production  of  the  product  as  a  function  of  the 
solute  concentration  and  the  total  number  of  radicals  present  when  steady 
state  is  reached.   Since  {[R]}  is  not  a  function  of  time  at  steady  state 
it  follov;s  that  the  total  rate  for  a  system  containing  n  spurs  will  simply 
be  n  times  the  rate  for  one  spur.   Letting  I  be  the  dose  rate  (energy 
absorbed  by  the  system  per  unit  mass  per  unit  time) ,  2w  be  the  energy 
required  to  produce  a  radical  pair,  y  be  the  density  of  the  system,  and 
a  be  the  ratio  of  the  number  of  unrecombined  radicals  in  a  spur  at  steady 
state  (steady  state  here  referring  to  the  time  at  which  the  radical 
recombination  reaction  becomes  negligible  in  comparison  to  the  other 
reactions)  to  the  number  initially  present,  we  have  for  the  total  rate  of 
production  of  product  species  per  unit  volume 


d  .    .  T       ayl,  , 
Ttl[?h     -\,-^S]°    •  (84) 
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where     t  =  time  since  beginning  of  irradiation. 

This  shows  that  for  low  dose  rates  the  rate  of  product  formation  is  propor- 
tional to  the  total  dose,  It,  received  up  to  time  t,s:hce  we  have  assumed  that  the  chain 
reactions  in  each  spur  are  unterminating.   In  actuality,  the  chains  will 
be  terminated  at  some  time,  giving  a  total  number  of  product  molecules 
for  each  spur.   For  this  case,  the  product  formation  rate  would  be 
proportional  to  the  dose  rate. 

Another  method  for  treating  this  case  which  is  somewhat  more  valid 
is  to  consider  the  reaction  as  consisting  of  a  diffusion-kinetics  portion 
followed  by  homogeneous  kinetics.   For  the  diffusion-kinetics  portion, 
if  t  is  the  time  at  which  ND(t)  levels  off,  at  which  time  we  have  the 
defined  quantities  Np  £  Np(t),  NR  =   NR(t) ,  NRg  =   N   (t)  .   if  e  is 
defined  as  the  energy  deposited  per  spur  then  the  number  of  spurs  per  cm 
per  second  is  equal  to  ul/e,  and  if  N.  is  Avagadro's  number,   so  the 
product  rate  of  formation  for  the  spur  reaction  is  given  by 


_l_i  =  103^__  Np  £  KNT  .  (85) 

dt       £NA        p 


For  the  subsequent  homogeneous  reaction  we  have  the  rate  equations 


d[R]  ,2 


dt 


=  KN     -  k      [&]z   -  k      [R][S]  +  k      [R    ][RS]    ,  (86) 

R  11  m  23      2 
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d[RS] 
dt 


-i  =  KNRC  +  k   [R] [S]  -  k   [R  ] [RS]  . 
Kb     ml  Jl  J     2  3   2 


(87) 


Since  at  steady  state  f*l!L  =  0  and  dfRS]  =  0  Eqs .  (86)  and  (87) 

dt  dt 

lead  to 


[R]  - 


—  <NR  +  NRS> 

c 

.  11 


-h 


(88) 


Substituting  this  into  the  product  rate  equation  given  analogously 
to  Eq.  (83)  but  for  homogeneous  reaction  as 


agL-k^SHt]  , 


(89) 


we  have 


d[P] 


dt 


■ku[S] 


(N   +  N   ) 

R    RS 


-11 


(90) 


The  overall  rate  is  then  equal  to  the  sum  of  the  rates  from  the 
homogeneous  reaction  and  the  spur  reactions 


dC!L.^T  +  k  [s] 

dt      P    ^ 


—(N  +  N   ) 
k    R    RS 

11 


(91) 


The  important  observation  from  Eq .  (91)  is  that  at  low  dose  rates, 
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the  product  rate  of  formation  is  equal  to  the  sum  of  a  term  linearly- 
dependent  on  the  dose  rate  and  a  term  with  a  square  root  dependence  on 
the  dose  rate  and  a  linear  dependence  on  the  solute  concentration.   This 

derivation  is  only  valid,  however,  where  the  depletion  of  R  and  S  is 

2 

negligible. 

For  the  case  of  high  dose  rates,  the  radiation  can  be  considered  as 
distributed  uniformly  so  that  homogeneous  kinetics  alone  apply.   The 
equations  describing  the  system  are  the  same  as  those  given  previously 
without  the  brackets.  When  steady  state  is  reached  in  this  case,  the 
time  rate  of  change  of  radical  concentration  and  the  time  rate  of  change 
of  intermediate  concentration  will  be  zero.   The  losses  in  radicals  will 
be  balanced  by  the  production  from  radiation  initiation  with  a  rate  ul/w. 
Setting  the  rate  of  radical  production  equal  to  zero 


ul/w  -  k^lR]2  -  k  [R][S]  +  k23[R2][RS]  =  0.  (92) 


Setting  the  intermediate  product  rate  of  formation  equal  to  zero 


kllf[R][S]  -  k23[R2][RS]  =  0.  (93) 


Combining  Eqs.  (80),  (92)  and  (93)  gives  (with  braces  deleted) 


d[P] 
dt     ^ 


yl 

wk 
11 


h 


(94) 
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In  this  case  we  see  that  the  rate  of  product  formation  is  proportional  to 
the  square  root  of  the  radiation  intensity,  and  is  also  proportional  to 
the  solute  concentration. 
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3.0  DISCUSSION  AND  RESULTS 

The  parameters  that  may  be  varied  in  the  solution  of  Eqs.  (65) 
through  (69)  are 

1)  geometry, 

2)  initial  solute  concentrations, 

3)  initial  free  radical  spatial  distribution, 

4)  rate  constants, 

5)  diffusion  coefficients, 

6)  initial  number  of  radicals. 

In  this  report  only  variations  of  initial  solute  concentrations,  rate 
constants,  and  initial  number  of  radicals  are  considered.   The  geometry 
is  taken  to  be  spherically-symmetric,  the  initial  free  radical  distri- 
bution is  taken  as  Gaussian  of  a  fixed  size,  and  the  same  value  is  used 
for  all  the  diffusion  coefficients.   Even  with  these  simplifications, 
the  variations  of  the  other  parameters  must  be  severely  restricted 
because  of  the  magnitude  of  computer  time  necessary  for  the  numerical 
solution  of  Eqs.  (65)  through  (69).   For  this  reason,  care  must  be  used 
in  selection  of  values  for  the  parameters  so  that  the  maximum  amount  of 
information  can  be  obtained.   A  table  of  the  parameters  used  in  this 
study  is  given  in  Table  I.   For  all  the  runs,  the  values  of  the  diffusion 
coefficients  were  taken  as 

D  =  2  x  10~5  cm2/sec,  (95) 
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and  the  value  of  the  parameter  representing  the  size  of  the  spur 
initially  is  taken  as 

R  =  7.07  X  .  (96) 

o 

These  particular  values  are  chosen  as  being  representative  of  those 

encountered  in  many  reaction  schemes. 

For  the  set  of  parameters  listed  in  Table  I,  a  digital  computer 
TABLE  I.   LIST  OF  RUNS  AND  THE  VALUES  OF  THE  PARAMETERS 


Run  Number 

k 
11 

k 

I  it 

k 
23 

N 
o 

[R   ] 
2    O 

[s] 

o 

1 

io-u 

io-  » 

10" 

11 

6 

IO"2 

IO"2 

2 

II 

It 

•  1 

ii 

10" l 

IO"1 

3 

M 

II 

II 

H 

1 

1 

4 

II 

II 

II 

ii 

10 

10 

5 

•  1 

II 

•  1 

ii 

io" " 

1 

6 

•  1 

II 

II 

n 

1 

10" ' 

7 

II 

II 

II 

9 

io' l 

II 

8 

II 

II 

II 

12 

II 

II 

9 

10" l3 

II 

II 

6 

II 

II 

10 

io-u 

io"  » 

II 

n 

•1 

II 

11 

II 

io'11 

10" 

9 

ii 

II 

II 

12 

II 

io"9 

10" 

11 

ii 

II 

II 

Note:   The  units  for  k  are  cm  /molecule-sec;  the  units  for  N  are 
radicals;  the  units  for  R2  and  S  are  moles/liter. 
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solution  of  Eqs.  (65)  through  (69)  lists  spatially-distributed  values  of 
the  concentrations  of  the  various  species  as  well  as  values  for  the  reac- 
tion losses  and  total  number  of  radicals  at  various  times.  As  the 
computer  output  contains  thousands  of  values,  the  best  method  of  represen- 
ting the  data  is  graphically.   Figures  1  through  12  show  plots  of  the 
following  quantities  versus  time: 

T 
2N   (t) :   twice  the  total  number  of  R  molecules  that  have  been 

R2 

formed  by  radical  recombination  up  to  time  t, 

T 

N   (t) :   the  total  number  of  RS  molecules  that  have  been  formed 
RS 

by  reaction  up  to  time  t, 

T 
N  (t) :    the  total  number  of  P  molecules  that  have  been  formed 

by  reaction  up  to  time  t, 

N  (t):    the  total  number  of  radicals  that  exist  at  time  t, 
R 

Z    (t) :    the  sum  of  the  reaction  loss  terms  and  the  number  of 
R 

radicals  existing  at  time  t  . 
Figures  13  through  16  show  sample  plots  of  the  spatial  variation  of 
the  concentrations  of  the  diffusing  species  as  a  function  of  time.   Note 
that  these  quantities  are  plotted  versus  the  dimensionless  quantity  p, 
rather  than  radial  distance.   Also  note  that,  since  all  of  the  diffusion 
coefficients  are  the  same  and  there  is  only  one  radical  species,  there  is 
a  one-to-one  correspondence  between  the  radial  position  determined  by  the 
same  value  for  P  for  each  plot. 

The  information  given  by  Figures  13  through  16  is  useful  in  under- 
standing the  physical  process  of  the  reaction  system.   Figure  13  shows 
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that  the  radicals  have  an  initial  Gaussian  distribution  which  is  maintained 
to  some  degree  with  a  decrease  in  amplitude  to  zero  as  the  radicals  diffuse 
and  are  lost  through  reaction.   During  this  time,  the  size  of  the  spur  is 
increasing  (remember  that  a  constant  value  for  p  corresponds  to  tiir.e- 
increasing  radial  distance)  to  many  times  its  original  size. 

Figure  14  shows  that  the  concentration  of  the  intermediate  product 
is  initially  zero,  builds  up  to  some  maximum  value,  and  then  returns  to 
zero.   The  buildup  is  caused  by  the  reaction  of  radicals  with  the  solute, 
S,  and  the  subsequent  decay  is  caused  by  diffusion  and  reaction  of  the 
intermediate  product  with  the  solute,  R  .   It  is  interesting  to  note  that 
the  shape  of  the  curve  is  approximately  Gaussian. 

Figure  15  shows  that  the  concentration  of  the  solute,  R  ,  is 
initially  constant  at  some  value,  builds  up  to  some  maximum  value,  and 
then  returns  to  its  initial  value.   The  buildup  is  caused  by  radical 
recombination,  and  the  decay  is  caused  by  the  reaction  of  R  with  the 
intermediate  product  and  by  diffusion.   This  spatial  distribution  is 
approximately  that  of  a  Gaussian  plus  a  constant. 

Figure  16  shows  that  the  concentration  of  the  solute,  S,  is 
initially  constant  at  some  value,  decreases  to  some  minimum  value  in  the 
spur  center,  and  then  returns  to  its  initial  value.   The  decrease  is 
caused  by  the  reaction  of  the  solute  with  the  radicals,  and  the  buildup  is 
caused  by  diffusion.   This  curve  has  the  approximate  shape  of  a  constant 
minus  a  Gaussian. 

Referring  to  Figures  1  through  12  it  is  possible  to  describe  the 


62 


results  in  a  general  fashion  by  discussion  of  the  following  eight  types 
of  parameter  variations 

1)  Increasing  both  solute  concentrations, 

2)  Increasing  the  concentration  of  R2, 

3)  Increasing  the  concentration  of  S, 

4)  Increasing  the  initial  number  of  radicals, 

5)  Decreasing  the  rate  constant  for  the  radical  recombination 
reaction, 

6)  Decreasing  the  rate  constant  for  the  radical-solute  reaction, 

7)  Increasing  the  rate  constant  for  the  chain  propagation 
reaction, 

8)  Increasing  the  rate  constant  for  the  radical-solute  reaction. 
The  effect  of  increasing  both  solute  concentrations  may  be  seen  from 

Figures  1  through  A.   The  effect  on  the  current  amount  of  radicals  present 
is  to  cause  a  more  rapid  initial  decrease  due  to  higher  reaction  rates 
with  the  solute,  S,  and  a  leveling  off  at  greater  amounts  of  radicals. 
The  effect  on  the  total  amounts  of  intermediate  and  final  product  formed 
is  to  cause  a  more  rapid  increase  due  to  the  higher  concentrations  of 
intermediate  product  achieved  coupled  with  the  higher  concentration  of 
solute,  R2.   The  effect  on  the  amount  of  R~  formed  by  radical  recombination 
is  to  cause  a  leveling  off  at  lower  amounts  due  to  the  decrease  in  the 
radical  concentration  caused  by  radical-solute  reaction. 

The  effect  of  increasing  the  concentration  of  solute  R  can  be  seen 
by  comparison  of  Figure  6  with  Figure  2.   The  effect  on  the  current  amount 
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of  radicals  present  is  to  cause  an  increase  due  to  the  higher  reaction 
rate  between  R  and  RS  causing  fewer  radicals  to  be  bound  up  with  the 
solute,  S.   The  effect  on  the  intermediate  and  final  product  is  to  cause 
an  increase  in  their  rates  of  formation  and  to  cause  the  difference 
between  their  respective  amounts  to  be  smaller,  this  latter  effect  being 
caused  by  a  more  rapid  reaction  rate  between  R2  and  RS.   The  effect  on  the 
amount  of  R~  formed  is  to  cause  an  increase.   This  is  the  result  of  a 
higher  concentration  of  radicals  present  since  fewer  are  retained  in  the 
form  of  RS  molecules. 

The  effect  of  increasing  the  concentration  of  solute  S  can  be  seen 
by  comparison  of  Figure  5  with  Figure  2.   The  effect  on  the  current  amount 
of  radicals  present  is  to  cause  a  more  rapid  decrease  to  a  lower  value 
before  leveling-off  is  achieved.   This  is  due  to  the  higher  reaction  rate 
between  the  radical  and  the  solute,  S.   The  effect  on  the  intermediate 
product  is  to  cause  a  large  increase  in  amount  formed,  the  same  being 
true  for  the  final  product  to  a  lesser  degree  since  the  amount  of  final 
product  has  a  higher-order  dependence  on  the  concentration  of  the  solute, 
S.   The  effect  on  the  amount  of  R  formed  is  to  cause  a  decrease.   This 
is  due  to  the  concentration  of  radicals  being  lowered  by  reaction  with  the 
solute,  S. 

The  effect  of  increasing  the  initial  number  of  radicals  can  be  seen 
by  comparison  of  Figures  2,  7,  and  8.   The  effect  on  the  current  amount 
of  radicals  present  Is  to  cause  a  similar  rate  of  decrease  to  a  relatively 
smaller  lcveling-off  value  (in  comparison  to  their  initial  values).   This 
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is  due  primarily  to  an  increase  in  the  radical  recombination  rate.   The 
effect  on  the  intermediate  and  final  product  is  small  due  to  the  high 
radical  recombination  rate  decreasing  the  amount  of  radicals  available 
for  chain  propagation.   The  effect  on  the  amount  of  R  formed  is  to  cause 
a  relative  increase  due  to  the  higher  radical  recombination  rate. 

The  effect  of  decreasing  the  rate  constant  for  the  radical 
recombination  reaction  can  be  seen  by  comparison  of  Figures  9  and  2. 
These  curves  show  the  expected  results  of  higher  amount  of  radicals  present, 
smaller  amount  of  IU  formed,  and  greater  amounts  of  intermediate  and  final 
product  formed,  this  last  being  the  consequence  of  higher  radical 
concentrations . 

The  effect  of  decreasing  the  rate  constant  for  the  radical-solute 
reaction  can  be  seen  from  Figures  10  and  2.   The  effect  on  the  current 
amount  of  radicals  is  to  cause  an  increase  since  the  radical-solute 
reaction  rate  is  slower.   The  effect  on  the  intermediate  and  final 
product  is  to  lower  their  values  and  cause  the  relative  differences  between 
these  products  to  be  smaller.   The  latter  is  caused  by  the  R-S  reaction 
being  slow  in  comparison  to  the  R?-RS  reaction.   The  effect  on  the  amount 
of  R~  formed  is  to  cause  an  increase  since  the  radical  concentration  is 
high  due  to  little  reaction  with  the  solute,  S. 

The  effect  of  increasing  the  rate  constant  for  the  chain  propagation 
step  can  be  seen  by  comparison  of  Figures  11  and  2.   The  effect  on  the 
current  amount  of  radicals  present  is  to  cause  an  increase  since  fewer  are 
being  retained  in  the  form  of  the  intermediate  product.  This  is  offset 
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somewhat  by  the  higher  concentrations  of  radicals  causing  more  radical 
recombination.   The  effect  on  the  intermediate  and  final  product  is  to 
increase  the  amounts  formed.   This  is  due  to  the  higher  radical  concen- 
tration causing  a  faster  chain  reaction.   The  major  effect  is  to  cause  the 
amount  of  P  formed  to  be  very  nearly  that  of  RS  formed  since  the  R--RS 
reaction  is  much  faster  than  the  R-S  reaction. 

The  effect  of  increasing  the  rate  constant  for  the  radical-solute 
reaction  can  be  seen  from  Figures  12  and  2.   These  curves  show  the 
expected  results  of  smaller  amount  of  radicals  present  and  R_  formed,  and 
greater  amounts  of  intermediate  and  final  products  formed. 

Having  formed  the  above  general  observations,  it  is  important  to 
obtain  a  qualitative  understanding  of  the  effects  that  spur  overlap  should 
have  on  the  rate  of  product  formation.   Here  we  will  only  consider  the 
case  where  overlap  occurs  before  the  amount  of  radicals  reaches  the 
leveling-off  region  shown  in  Figures  1  through  12,  as  little  radical 
recombination  occurs  after  this  time.   Also  necessary  is  a  knowledge  of 
the  spur  size  increase  with  time  so  that  one  can  estimate  the  spur 
separation  distances  that  would  be  necessary  for  the  effects  of  overlap 
to  become  important.   Table  II  gives  approximate  values  for  the  factor 
by  which  the  Gaussian  width  at  half-maximum  has  increased  as  a  function 
of  time.   These  values  were  determined  from  the  computer  solution  data 
for  Run  1  and  are  fairly  representative  values  for  all  the  runs.   The 
deviations  from  these  values  for  the  other  runs  does  not  exceed  approx- 
imately five  percent. 
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TABLE  II.   SPUR  SIZE  AS  A  FUNCTION  OF  TIME 


Time 

(seconds) 

r^(t)/rii(o) 

1.5  x 

io-n 

1.1 

1.8  x 

10"10 

1.7 

1.4  x 

10"9 

3.6 

1.2  x 

10"8 

10 

1.1    X 

ID"7 

30 

The  basic  effects  that  overlap  will  have  are  as  follows.   The  over- 
lapping of  the  radicals  in  the  spurs  will  cause  the  recombination  reaction 
to  increase,  decreasing  the  radical  concentration  and  increasing  the  R? 
concentration.   The  decrease  of  the  radical  concentration  will  cause  a 
decrease  in  the  rate  of  production  of  the  intermediate  product.   It  would 
be  expected  that  initially,  the  rate  of  production  of  the  final  product 
would  increase  because  of  the  increased  concentration  of  the  intermediate 
product  in  the  overlap  region,  along  with  the  increased  concentration  of 
R_  from  the  increased  recombination  occurring.   The  rate  will  subsequently 
fall  below  the  value  without  overlap  because  of  depletion  of  available 
intermediate  product. 

Having  established  that  the  basic  effects  of  spur  overlap  depend 
primarily  on  the  radical  recombination,  the  important  question  is  then 
how  the  presence  of  a  chain  reaction  affects  the  concentration  of  the 
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radicals  as  a  function  of  position  and  time.   In  order  to  simplify  the 
analysis  it  is  necessary  to  assume  that  the  spatial  concentration  profile 
of  the  radicals  remains  of  a  Gaussian  shape  at  all  times  so  that  the 
important  information  is  represented  by  values  of  the  radical  concen- 
tration at  the  half-maximum  position.   This  assumption  is  quite  good  for 
the  runs  made  in  this  report.   The  values  obtained  are  listed  in  Table  III 
and  plotted  in  Figure  17.   The  concentrations  have  been  normalized  to 
their  respective  initial  values  in  the  center  of  the  spur. 

In  order  to  obtain  an  upper  limit  for  these  curves  for  the  case  of 
no  radical  recombination  and  infinite  propagation  rate  (that  is,  for  k^ 
■  0  and  k23  =  m)   we  note  that  this  would  correspond  to  the  case  for  radical 
diffusion  with  no  reaction  for  which  we  have  the  analytical  solution 


[R(r,t)]  = 


No 


3/2 


P' 


p(2lT+4Dt)  I 


2R2+ADt 
o 


(97) 


Noting  that 


[R(o,o)l  =  [R]   = 


No 


[2< 


3/2 


(98) 


wc  have  for  Eq.  (97) 


[R(r?jtt)1 


2R6* 


2R2+4Dt 
o 


3/2  -r2/  [2R2+4Dt 
e   *    o 


(99) 
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TABLE   III.      RADICAL   CONCENTRATIONS   AT   HALF-MAXIMUM  POSITION   AS 
A  FUNCTION   OF  HALF-MAXIMUM  POSITION 


RUN 

[R],/[R] 

*2                O 

Vt) 

1 

0.376     0.796*10"1     0.619*10""2     0.175X10"3     0.480X10"5 

9.21       14.4                  30.0                  81.4                  246 

2 

0.374     0.751*10-1      o.462x10-2     0.158xl°~3     0.544x10~5 

9.20        14.3                  29.2                   81.0                  248 

3 

0.350     0.564x10~l     0.451x10"2     0.158x!0~3 

9.20       13.6                  29.0                  90.3 

4 

0.230     0.566x10-!     0.812x10-2 

9.05        13.5                  23.9 

5 

0.348     0.399x10-!     0.814x10"3     0.348xl°_lt 

9.14        13.3                  29.6                   82.6 

6 

0.374     0.778x10-!     0.634xl°-2     0.254x10~3 

9.21        14.5                  29.8                  81.7 

7 

0.352     O.640x10-1     0.399x10"2     0.136x10~3     0.460x10~5 

9.40        14.7                  29.4                  81.1                  251 

8 

0.334     0.557x10-!     0.354x10"2     0.120x10-3     0.386x10~5 

9.56        15.1                  29.6                   81.2                  260 

9 

0.424     0.120                0.774x10~2     0.245xl°~3     0.358x10"5 

8.84       13.0                  28.0                  80.4                  254 

10 

0.376     0.801x10-!     0.652x10~2     0.262x10~3     0.870x10~5 

9.21        14.5                   30.1                  82.1                  256 

11 

0.375     0.801x10-!     0.653x10-2 

9.20        14.5                   30.1 

12 

0.346     0.682X10-1 

8.75       11.7 
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where     r^  =  radial  position  at  a  concentration  equal  to  half  the 

maximum  value. 
To  obtain  ri  (t)  we  have 


[R(r,  ,t)l 
3 

[R(o,t)l 


iL 


=  h   -  e 


2R£+4Dt 


(100) 


so  that 


r? (t)  =  [2R*  +  4Dt]ln(2)  . 


(101) 


Substituting  Eq.  (101)  into  Eq.  (99)  gives 


[R] 

[Rio 


=  0.816 


T=Xh 


if* 


(102) 


[Rl. 

Equation  (102)  is  the  relation  for  the  maximum  1  that  could  be 

[R]0 
obtained  from  the  chain  reaction  scheme  considered,  for  any  value  of 

ri  (t) .   This  equation  is  plotted  as  the  dashed  line  in  Figure  17.   A 

log-log  plot  was  used  in  order  to  investigate  deviations  of  the  numerical 

solution  values  from  a  straight  line.   The  only  serious  deviations 

occurred  for  Runs  3  and  A  so  these  were  omitted  from  the  plot.   The 

deviations  from  a  straight  line  are  caused  by  deviations  in  the  spatial 

distribution  of  the  concentration  of  radicals  from  a  Gaussian  form.   The 
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form  was  flattened  in  the  center  of  the  spur  for  these  cases  due  to  the 
high  reaction  rates  caused  by  large  solute  concentrations. 

In  order  to  draw  useful  conclusions  from  Figure  17  we  must  keep  in 
mind  that  the  purpose  here  is  to  attempt  to  determine  qualitatively  whether 
reaction  can  increase  the  importance  of  overlap  effects.   Thus,  the 
larger  [R]i7[R]   is  for  a  given  n  ,    the  more  important  overlap  effects 
become.  We  notice  from  Figure  17  that  the  run  which  would  have  least 
importance  for  overlapping  is  Run  12.   Comparing  Run  12  with  Run  2,  and 
comparing  the  parameters  that  produced  these  runs,  we  notice  that  we  could 
interpret  Run  2  as  the  chain  reaction  and  Run  12  as  the  non-chain  reaction. 
The  reason  why  it  is  possible  to  treat  one  reaction  as  a  chain  reaction 
and  the  other  as  the  non-chain  reaction  is  basically  that  in  one  reaction, 
the  chain  propagation  step  is  more  pronounced  than  in  the  other  reaction. 
Thus,  although  both  reactions  are  actually  chain  reactions,  it  is  feasible 
for  comparative  purposes  to  term  the  one  with  the  greater  chain  length 
the  chain  reaction  and  the  other  the  non-chain  reaction.   Of  course,  the 
actual  non-chain  reaction  would  be  the  one  for  which  the  reaction  rate 
constant  for  the  chain  propation  step  was  identically  zero,  but  this 
would  reduce  the  problem  to  one  of  recombination  only.   Since  the  values 
for  Run  2  from  Figure  17  are  higher  than  those  for  Run  12,  the  conclusion 
is  that  overlap  effects  should  be  more  important  for  chain  reactions  than 
for  non-chain  reactions  due  to  the  higher  concentrations  of  radicals 
present.   A  more  realistic  comparison  can  be  obtained  by  interpreting 
Run  2  as  the  non-chain  reaction  and  Run  11  as  the  chain  reaction.   This 
comparison,  also,  shows  that  the  radical  concentrations  are  higher  for  the 
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chain  reaction  case. 

The  runs  which  gave  the  highest  radical  concentrations  are  1,  6,  9, 
10,  and  11.   It  is  not  possible  to  draw  general  conclusions  from  these 
runs  except  to  note  that  in  comparison  to  Run  2,  higher  concentrations 
were  obtained  by  decreasing  the  recombination  rate  constant;  by 
decreasing  the  radical-solute  rate  constant;  and  by  increasing  the  chain 
propagation  rate  constant. 

An  interesting  feature  is  noticed  by  comparing  Runs  2,  7,  and  8. 
These  cases  correspond  to  an  increase  in  the  initial  number  of  radicals 
in  the  spur  and  it  can  be  seen  that  although  there  was  some  decrease  in 
relative  concentration  when  the  initial  number  of  radicals  was  increased 
from  6  to  9,  there  was  no  further  decrease  in  relative  concentration  when 
the  number  was  increased  from  9  to  12  radicals  per  spur.   Therefore  it  may 
be  concluded  that  the  effect  of  increasing  the  radical  concentration  has 
only  a  minor  effect  on  the  relative  concentration  as  a  function  of  width 
at  half-maximum,  although,  of  course,  the  actual  radical  concentration  is 
higher. 

A  further  substantiation  that  a  chain  reaction  significantly  increases 
the  radical  concentration  is  seen  by  comparison  of  Run  5  with  Run  6. 
Because  of  the  initial  concentrations  involved,  Run  5  favors  the  solute 
reaction  step  and  Run  6  favors  the  chain  propagation  step,  so  that  we  may 
interpret  Run  5  as  the  non-chain  reaction  and  Run  6  as  the  chain  reaction. 
Once  again  we  see  that  effect  of  a  chain  reaction  is  to  significantly  in- 
crease the  radical  concentration,  so  that  overlap  will  be  more  important. 
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4.0  ACCURACY  OF  THE  FINITE  DIFFERENCE  SOLUTION 
AND  ANALYSIS  OF  ERRORS  FOR  THE  RUNS 

In  order  to  estimate  the  accuracy  of  "and  determine  the  sources  of 
error  in  the  finite  difference  solution  values  two  other  runs  were  per- 
formed.  The  first  was  a  solution  of  the  two-radical,  non-depletion  of 
solute,  model  used  by  Dyne  and  Kennedy  (A)  with  the  reaction  scheme  given 
as 


R  +  R  — >  R   Rate  Constant  =  It,,  ,  (105) 
2                 RR 

R  +  R*  — >  RR*   Rate  Constant  =  k^*  ,  (106) 

R*  +  R* — >  R*   Rate  Constant  =  kR*R*  ,                   (107) 


R  +  S — >  RS   Rate  Constant  =  kRS  ,  (108) 


R*  +  S >  R*S   Rate  Constant  =  k  (109) 

R*S  . 


The  values  of  the  parameters  used  in  the  solution  were 

D  =  «  , 
S 

k   -  10"11  cm3/sec-radical  , 
RS 

—  1 1    3 
k   ■  1.2  x  10    cm  /sec-radical  , 
RR 

k   .  =  3.0  x  10"11  cm3/sec-radical  , 
RR* 

k  .  .  =  0.9  x  10-11  cm3/sec-radical  , 
R*R*  * 

D  »=  8.0  x  10" 5  cm2 /sec  , 
K 
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DR*=  2.0  x  10~5  cm2/sec  , 

N  =  6.0  radicals  , 
o 

N*   6.0  radicals  , 
o 

o 

R  =  7.07  A  , 
o 

R*  -  7.07  A  , 
o 

[S]  =  10~3  mole/liter  . 
o 

The  results  of  the  computer  solution  compared  to  those  published  by 
Dyne  and  Kennedy  and  those  computed  by  F.  E.  Haskin   (6)  using  the  pre- 
scribed diffusion  hypothesis  are  listed  in  Table  V.   In  this  table,  N 

and  N*  are  the  initial  number  of  R  and  R*  radicals  respectively,  N   , 

°  R2 

NR*,  and  N  *  are  the  final  numbers  of  R2,  Ro*>  and  RR*  molecules  formed 

by  radical  recombination  respectively,  and  N  and  N   are  the  final 
J  R      R* 

numbers  of  radicals  which  reacted  with  the  solute. 


TABLE  V.   COMPARISON  OF  THE  TWO-RADICAL  MODEL  RESULTS 


Quantity 

This  Work 

Dyne 

and  Kennedy 

Haskin 

2N  /N 
R2   o 

0.115 

0.093 

0.127 

2NR*/NJ 

0.217 

0.183 

0.233 

NRR*/No 

0.389 

0.417 

0.431 

VNo 

0.496 

0.490 

0.439 

NR*/N* 

0.394 

0.400 

0.333 
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Since  Dyne  and  Kennedy's  results  are  supposed  to  be  accurate  to  1% 

or  so,  it  is  seen  that  for  the  number  of  spatial  mesh  points  used  in  this 

work,  only  an  accuracy  of  10%  or  so  can  be  expected.  A  higher  accuracy 

was  not  sought  after  in  this  work  since  it  was  not  the  numerical  value 

that  was  important  but  the  trend  of  the  solutions.   To  obtain  a  higher 

accuracy,  a  smaller  mesh  size,  and  therefore  more  mesh  points,  would  have 

to  be  used,  thus  significantly  increasing  computation  time.   For  instance, 

if  the  mesh  size  in  the  spatial  coordinates  were  halved,  the  computation 

time  would  have  increased  by  a  factor  of  eight  since  the  stability 

relation  forces  a  factor  of  four  decrease  in  time  mesh  size.   This  added 

computation  time  was  not  deemed  necessary  in  this  work.   Another  estimate 

of  the  accuracy  is  given  by  Kuppermann  (7)  as  being  the  percent  error  in 

the  £   calculated  for  the  runs.   Kuppermann  states  that  in  the  runs  he 

performed,  the  percent  error  in  £  was  a  good  estimate  of  the  percent 

K 

error  in  the  values  determined.   For  the  run  described  above,  the  percent 
error  in  the  E_  value  was  about  10%,  tending  to  agree  with  the  numerical 
results  listed  in  Table  V. 

In  order  to  show  that  these  magnitudes  of  errors  do  not  alter  the 
general  shapes  of  the  time  plots,  as  well  as  indicating  the  validity  of 
using  the  prescribed  diffusion  hypothesis  even  in  this  complex  reaction 
system,  the  results  obtained  from  this  work  are  plotted  in  Figure  18  along 
with  those  obtained  by  Haskin   (dashed  lines).   It  is  seen  that  even 
though  the  percent  error  in  the  final  values  is  large,  the  general  shapes 
of  the  curves  do  not  change  much  in  a  qualitative  sense,  and  it  is  these 
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trends  that  are  important  in  this  work. 

Since,  in  general,  the  only  method  for  locating  the  sources  of  errors 
in  any  numerical  solution  is  to  change  the  mesh  sizes,  Run  3  was  recal- 
culated in  Run  3-A  using  a  time  mesh  size  that  was  decreased  by  a  factor 
of  two.   Since  the  results  obtained  were  so  close  to  the  original,  it  is 
necessary  to  show  the  differences  in  tabular  form.   The  results  are  given 
in  Table  VI. 


TABLE  VI.   EFFECT  OF  THE  TIME  MESH  SIZE  ON  THE  NUMERICAL  RESULTS 


Problem  Time              Quantity                Run   3                  Run   3-A 
(sec) 

8.35  x  10"11                   NRS(t)                   2'17                         2'20 

Np(t)                0.522                      0.533 

2NR   (t)                0.775                      0.782 

NR(t)                   3.58                         3.58 

ER(t)                   6.00                         6.01 

5.55  x  10"10                  N5s<t)                  9.08                        9.04 

Np(t)                   6.73                         6.70 

2Nr    (t)                  1.18                        1.19 

NR(t)                   2.36                         2.35 

ER(t)                  5.89                        5.88 

5.35   x   10~9                     NRs(0                   70-°                         69-9 

Np(t)                   67.8                         67.6 

2N'J    (t)                   1.40                         1.40 

NR(t)                   2.11                         2.07 

ErCO                  5.77                        5.73 
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Since  this  table  shows  that  the  decrease  in  the  time  mesh  size  had  little 
effect  on  the  results,  the  source  of  error  must  lie  in  the  spatial 
portion  of  the  finite  difference  solution.   The  two  possible  regions  of 
error  here  are  the  effect  of  only  using  a  finite  portion  of  the  Gaussian 
and  the  finite  difference  error  in  the  computation  of  the  spatial  deri- 
vatives.  Since  the  program  used  includes  more  than  99.99%  of  the 
Gaussian  in  its  calculations,  the  main  source  of  error  must  lie  in  the 
computation  of  the  derivatives.   This  is  a  reasonable  conclusion  since 
only  thirty  spatial  mesh  points  were  used.   Therefore,  to  obtain  higher 
accuracy,  it  is  necessary  to  use  a  larger  number  of  spatial  mesh  points. 

Concerning  the  analysis  of  errors  in  the  runs  there  is  only  one 
major  correction  that  could  be  made,  and  this  is  for  the  curves  showing 
the  time  behavior  of  the  current  amounts  of  radicals  present  at  any 

time.   If  one  plots  I  and  N  (t)  for  any  run  on  an  expanded  scale  near 

R      R 

the  end  of  the  runs  as  is  done  in  Figure  19  for  Run  2,  it  is  noticed  that 

their  respective  slopes  are  very  nearly  the  same.   If  one  then  plots 

T        T 
N   (t)  and  N  (t)  for  any  run,  as  was  done  in  Figure  20  for  Run  1,  it  is 

Kb  A 

seen  that  these  quantities  are  very  linear  with  time  near  the  end  of  the 

computation.   Since  [S]  at  the  times  considered  is  quite  constant,  it  is 

seen  from  Eq.  (83)  that  in  order  for  the  time  rate  of  change  of  product 

formed  to  be  constant,  then  N  (t)  must  be  constant.   If  we  then  can 

R 

assume  that  N  (t)  is  nearly  constant,  we  have  the  conclusion  that  the 
R 

major  portion  of  the  slope  of  I  is  due  to  spurious  losses  of  radicals, 

R 

that  is,  losses  due  to  the  numerical  method  employed.  Therefore,  in  all 
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the  runs  made,  N  (t)  is  actually  much  more  constant  towards  the  end  of 
R 

the  runs  than  indicated  by  the  figures. 
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5.0   CONCLUSIONS 

The  results  of  this  work  have  led  to  the  major  qualitative  con- 
clusion that  effects  of  spur  overlap  will  be  important  at  lower  dose 
rates  for  chain  reactions  than  for  non-chain  reactions.   This  conclu- 
sion is  based  on  the  fact  that  for  the  runs  made  for  the  reaction 
scheme  considered,  there  was  a  significant  increase  in  the  radical 
concentration  throughout  the  history  of  the  spur  for  chain  reactions 
as  opposed  to  non-chain  reactions,  indicating  that  radical  interactions 
between  neighboring  spurs  would  be  greater  for  the  same  dose  rates  for 
the  chain  reaction.   It  has  also  been  shown  that  the  overall  reaction 
rate  for  the  system  considered  should  be  lower  than  that  predicted 
from  homogeneous  kinetics  for  the  case  of  high  dose  rates.   This 
decrease  in  the  overall  reaction  rate  is  caused  by  radical  recombin- 
ation in  each  spur  before  complete  spatial  uniformity  of  the  radical 
concentrations  is  reached.   This  inhomogeneity  will  also  cause  the 
overall  reaction  rate  to  depend  on  more  of  the  reaction  parameters 
than  predicted  from  homogeneous  kinetics. 

One  other  conclusion  that  was  reached  from  an  examination  of  the 
space-time  histories  of  each  of  the  species  for  each  of  the  runs  is 
that  a  modified  form  of  the  prescribed  diffusion  hypothesis  should 
give  good  results  for  the  reaction  scheme  considered  except  for  very 
high  solute  concentrations,  where  the  original  Gaussian  form  is 
flattened  in  the  center  of  the  spur. 
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6.0   SUGGESTIONS  FOR  FURTHER  STUDY 

This  study  was  intended  to  serve  as  a  basis  to  establish  the  impor- 
tance of  spur  overlap  for  chain  reactions  and  the  only  calculations 
performed  were  those  for  a  single  cluster  of  reacting  species.   This 
work  can  be  extended  by  obtaining  numerical  solutions  for  the  case  of 
spur  overlap  for  a  chain  of  equally-spaced  spurs  with  the  addition  of 
one  more  spatial  variable,  thus  creating  a  partial  differential  equation 
in  three  variables  to  be  solved.   It  is  estimated  that  increase  in 
computation  time  for  this  case  would  be  a  factor  of  ten  to  thirty  over 
that  for  the  two-variable  equation  used  in  this  work.   This  amount  of 
computation  time  is  not  prohibitive  if  the  parameters  used  in  the  study 
are  chosen  carefully. 

A  further  extension  of  this  work  is  to  perform  the  same  calculations 
used  in  this  work  with  the  use  of  a  modified  form  of  the  prescribed 
diffusion  hypothesis.   This  modification  would  involve  the  use  of  a 
Gaussian  shape  for  the  radicals  throughout  their  history,  a  constant  plus 
a  Gaussian  for  the  shape  of  the  species  formed  by  radical  recombination, 
a  constant  minus  a  Gaussian  for  the  solute,  and  a  Gaussian  for  the 
intermediate  product.   This  type  of  system  would  reduce  the  solution  of 
the  partial  differential  equation  to  the  solution  of  an  ordinary 
differential  equation,  amounting  to  a  tremendous  savings  in  computation 
time. 
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A.O  APPENDIX 


A.l  Justification  for  the  Variable  Transformation 


r 


n>4  (ti 


i(t0+t) 

Consider  the  case  where  free  radicals  diffuse  outward  from  a  point 
source  with  no  reaction.   The  diffusion  equation  is  given  as 


3C.  (r,t)     32C,(r,t)   2   3C.  (r,t) 

— i =  D  i +  ^d  — l , 

1     o      r  L 
3t  3r2  3t 


(110) 


with  boundary  conditions 


C,(r,o)  =  Ni6(r)/^Trr2,  (HI) 

i        o 


C^KO  =  0  ,  (112) 


3Mo,t) 

1   '   ■  0,  (113) 

3r 


where     6(r)  =  Dirac  delta  function, 
r  =  radial  distance, 
C.(r,t)  =  concentration  of  free  radicals  of  species  i, 
t  =  time, 
N  =  number  of  free  radicals  present  of  species  i. 
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The  solution  to  Eq.  (110)  with  the  given  boundary  conditions  yields 


NVr2/ADit 
Ct(r,t)  =- 


UttD^] 


3/2 


(114) 


If  the  time  scale  is  shifted  so  that  time  zero  in  a  new  scale  is 

equal  to  time  t   in  the  old  scale  and  a  new  quantity,  R  ,  is  defined  by 
o  o 


4D 


t*  =   2[R*]2   Eq.    (114)    becomes 


C^r.t)   = 


NW/UDit  +  2(rJ)21 


[7i{4Dit  +  2(R0)2}] 


3/2 


(115) 


which  has  the  initial  condition 


I  -r2/2(l62 


Ci(r,o)  = 


Noe 


[2^(rJ)? 


3/2 


(116) 


Eq.  (115)  therefore  gives  the  solution  to  the  diffusion-kinetics  equation 
in  the  absence  of  reaction  for  the  type  on  initial  condition  considered 
for  free  radicals  in  this  report. 

Defining  the  new  variable  N^(R)  as  the  number  of  radicals  present  in 
the  volume  bounded  by  r  =  R 


N^R) 


47ir2drCi(r,t) 


(117) 
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Defining  the  variable  transformation 


P 


r 


i   ; : (H8) 


/H<si  +  O 


Substituting  Eq.  (118)  and  Eq.  (115)  into  Eq.  (117)  gives 
N*    ?     „      N*p3 


"i<p>  -— %7i  UpiPi2=   -^-37i* 

(2»)  /2  o        24(tt)3/Z 


(119) 


The  result  obtained  in  Eq.  (119)  shows  that  a  fixed,  finite  range  in 
p  produces  a  time-independent  determination  of  the  number  of  radicals  for 
the  case  of  Gaussian  initial  conditions  and  absence  of  reaction.   There- 
fore the  introduction  of  this  transformation  into  the  diffusion-kinetics 
equation  will  allow  a  fixed,  finite  range  in  P  to  include  the  same 
amount  of  information  about  the  free  radical  concentration  at  every  time 
step.   It  can  also  be  shown  that  the  same  type  of  result  is  obtained  for 
the  case  of  cylindrical  tracks  with  the  transformation  given  by  Eq.  (118). 
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A. 2  Explanation  of  the  Computer  Program  Used  in  This  Work 


The  computer  program  discussed  here  numerically  solves  the  diffusion- 
kinetics  equations  for  spherically-symmetric  geometry  with  either  Gaussian- 
initial  free  radical  spurs  and  diffusing  solutes  or  Gaussian-initial  free 
radical  spurs  and  non-diffusing  solutes.   The  program  is  written  in  Fortran 
language  for  the  CDC-3600  computer  and  the  output  data  include   concentra- 
tion profiles  and  values  of  the  reaction  integrals  at  spaced  time  intervals 
during  the  calculation.   This  program  is  separated  into  several  subprograms 
for  ease  of  understanding.   The  following  list  indicates  the  program 
division: 

DIFFUSE  Main  Program 

SPURCYL 

INTMULT 

INTG1D 

FATES 


.  .  .  .  Subroutines 


DIFFUSE  Program 
The  function  of  this  program  is  to  read  and  write  the  input  data, 
initialize  all  variables  used,  calculate  initial  free  radical  and  solute 
concentration  profiles,  and  determine  the  spatial  and  time  increments 
such  that  stability  of  the  finite  difference  solution  is  initially  assured. 
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Input  Data 


Symbol 


Explanation 


NUMBER 

TTIME 

NUMRAD 
NUMSOL 
NUMRPTS 
ITYPESOL 

RZERO(I) 
D(I) 

NRXNTERM(I) 


IC(I,J), 
JC(I,J) 


S(I,L,M) 
C(1,I) 


Number  of  separate  solutions  to  the  diffusion-kinetics 
equations  to  be  run 

Length  of  time  considered  for  the  reaction 

Number  of  free  radical  species 

Number  of  solute  species 

Number  of  spatial  grid  points  used 

Code  number  for  type  of  solute  diffusion 

■  1  for  solute  diffusion 

=  0  for  lack  of  solute  depletion  assumption 

Initial  Gaussian  width  at  half-maximum  for  species  I 

Diffusion  constant  for  species  I 

Number  of  reaction  terms  appearing  in  the  diffusion- 
kinetics  equation  for  species  I 


Index  indentifiers  associated  with  the  reaction  terms. 
For  example,  if  in  the  diffusion-kinetics  equation  for 
species  1   the  fourth  reaction  term  involves  reaction 
between  species  2  and  species  3,  the  values  of  the 
index  indentifiers  become  IC(1,4)  =  2  and  JC(1,A)  =  3. 

Rate  constant  for  the  reaction  of  species  L  with  species 
M  in  the  diffusion-kinetics  equation  for  the  I   species 

Initial  number  of  free  radicals  of  species  I  present  or 
initial  concentration  of  the  Itn  solute  in  atoms/unit 
volume 


93 


Other  Important  Symbols  Used 


Symbol 


Explanation 


C(J,I)      Spatial  concentration  profiles  for  species  I  at 
position  J 

DELTAR      Spatial  increment,  Ap 

DELTAT      Time  increment,  Ax 

WNR(J)      Weight  factor  to  be  used  in  evaluation  of  volume 
integrals 

TAU(I)      Initial  time  parameter,  t1 


The  subprograms  used  by  this  program  are: 

SUBROUTINE  FATES  (1,  NUMRPTS ,  RADIAL,  WNR) 
SUBROUTINE  SPURCYL 
These  subprograms  are  explained  in  a  later  section. 


Logic  Diagram  for  Program  DIFFUSE 
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f      START  \ 


READ 
INPUT 
DATA 


Generate 
C(J,I) 

for  RADICALS 


WRITE 
input 
data 


Generate 
DELTAT 


Generate 
TAU(I) 


Generate 
DELTAR 


RADTAL  (J) 

«  JADELTAR 


-c 


CALL 
FATES 


for  Solute 

«=  .-1 


Adjust  DELTAT 
for  REACTION  RATE 
stability 


Generate 
C(J,I) 

for   SOLUTES 

G 


CALL 
SPURCYL 


> 
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no 


DECREASE 
NUMBER 
by  One 


STOP 
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DIMENSION 

DIMENSION 

DIMENSION 

DIMENSION 

COMMON  CZL: 

COMMON 

COMMON 

COMMON 

COMMON 


PROGRAM  DIFFUSE 

DIMENSION  RZERO(5) »D( 5  J  tC( 50  •  10    )  »NRXNTERM(5) »IC(5»5),JC(5,5) 
S(5»5»5)»A(5),  CZERO( 5)  »CLEFT(5) ,RADIAL( 50) 

WNR(50) »WTR(5»50 )  »AXIAL<  50)  ,WNZ(50)  »MARG<  5 »5 ) »CZLFFT ( 5 ) 
C0EFF1( 5)»C0EFF2(5)  ,C0EFF3( 5)  »C0EFF4 ( 5  0 , 5 ) »C0EFF5 (5*50) 
SUM(5) »CMULT(5»50) ,CZT|_OSS  (  5  »  5  ) »CLOSS ( 5 ♦ 5 ) »DLTM ( 5 ) 
SS(5,5) 
TAU(5 ) 
TIMFB( 5) 
TTIME 

SPUR, NUMR AD »NUMSOL»NUMRPTS,RZERO,D»NRADSOL,C»NRXN TERM, IC.JC 
ltS»NUMZPTS,SPURSEP» ALPHA, TIN  IT, DEL TAZ,A»DELTAR,DELTAT,CZERO,C LEFT , 
1Q,RADIAL,WNR,WTR,AXIAL,WNZ,PI , I TYPESOL ,MARG ,COEFF 1 , C0EFF2 ,C0EFF3 , 
1C0EFF4,C0EFF5-,TIME,DELTIME,SUM,CMULT,CZTL0SS,CL0SS,DLTM,CZLEFT 

1  FORMATUE15.8) 

2  FORMAT ( 1CI5) 

3  FOF  MAT( I5»E15.8) 

4  F0RMAT(2X,8E16.8) 
READ    2, NUMBER 
PRINT    2, NUMBER 

DO    643    1  =  1  ,10 
DO    643    J=l,50 

643  C( J, I )  =  0. 

DO    644    1=1 ,5 
DO    644    J=l ,5 

644  CLOSS( I ,J)=0. 
READ    1,TTIM1 

PRINT    1,TTIME  • 

NUMZPTS=1 

RFAD  2*NUMRAD,NUMS0L,NUMRPTS 

PRINT  2,NUMRAD,NUMS0L,NUMRPTS 

READ  2»ITYPES0L 

PRINT  2.ITYPES0L 

READ  1,  (RZERO( I ) , 1  =  1 ,NUMRAD) 

PRINT  1»(RZER0( I )  ,  I  =  1,NUMRAD) 

NRADSOL=NUMRAD+NUMSOL* I TYPESOL 

READ  1,  (D(  I  )  ,1  =  1  ,NRADSOL) 

PRINT  1, ( D( I ) , I=1,NRADS0L) 

READ  2*  (NRXNTERM(  I )  ,  I  =  1  ,NRADSOL ) 

PRINT  2, (NRXNTERM(  I )  ♦  I  =1  ,NRADSOL ) 

DO  100  I=1,NRADS0L 

N=NRXNTERM{ I ) 

DO  100  J=1»N 

READ  2»  IC(  I  ,J) ,JC( I ,J) 

PRINT  2,  IC( I  ,J) ,JC(  I ,J) 

II  =  IC(  hJ) 

JJ=JC( I *J) 

READ  1,  S( I ,1  I ,JJ) 

PRINT  4,S(  I , I  I  , JJ) 

NRADSOL=NUMRAD+NUMSOL 

READ  1 » (C(  1* I   )  » 1  =  1  »NRADSOL) 


100 
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PRINT    1  » ( C (  1  1 1        ) » I  =  1*NRADSCL) 

ALPHA=2. 

Z=NUVRPTS-1 

DEI  TAR=( ( 55.  )**0.5)/Z 

DEUAT  =  0.2*(DELTAR**2 ) 

DC  112  I=1»NUMRAD 

CZERC( I  )=C(  1  ♦  I   ) 

112  CLEFT( I )=CZERC( I ) 
DC  108  J=1.NUMRPTS 
Z  =  J-1 

108  RADIALf J) =Z*DELTAR 

CALL  FATES ( 1 »NUMRPTS t RADIAL »WNR ) 

PI=3. 1415926 

NR  =  NUf<'RPTS-l 

DC  113  I=1»NUMRAD 

C(NUMRPTS» I )=0. 

TAU( I )=(RZFRC(  I  )**2  )  /(2.*D(  I  )  ) 

DC  113  J=1,NR 

Z=J-1 

113  CCJtl   )=(CZERC( I  )/(  (2.*PI*(RZERC(  I  )**2. ) )**1.5) )*EXPF(  ((Z* 
1DELTAR)**2  »/4«  ) 

N=NUMRAD+1 

IF(NUMSCL.EQ.O)  GO  TC  177 

NBC=NUMRPTS+2 

DC  114  I=N.NRADSCL 

TAU( I )=TAU( 1 ) 

IE=I+NRADSCL 

DC  114  J=2*NBC 

C(JtIE)=C( 1» I ) 

114  C(JtI)«CCltI) 

DC  106  I=1»NRADSCL 

IF(  I.GT.NUMRAD.AND. ITYPESCL.EQ.O)  GC  TC  180 

N=NRXNTERM( I ) 

DC  106  L=l  iN 

I  I  =  IC( I »L) 

JJ^JC( I »L) 

IF(S( I ♦ I  I , JJ) .GE.O. )  GC  TC  106 

IF(Cdtl)  .FO.O.  )  GC  TC  106 

IF(C( 1 » I  I  ) .EO.O. )  GC  TC  106 

!FtCfl«JJ) .EO.O. )  GC  TC  106 

G  =  0.i*c(  1 »  I  ) /(TALK  1 )*A8SF(S(  I  ♦ I  I  »JJ)  )  *C  (  1  .  I  I  )*C( ltJJj  ) 

IF(DELTAT.GT.G)  DELTAT=G 
106  CONTINUE 
180  CONTINUE 
177  CONTINUE 

PRINT  4,DELTAZ»DELTAT,T INI T 

PRINT  4iD£LTAR 

PRINT  4»((C(J.I   )  iJ  =  l.NUMRPTS)  .1  -ltNRADSCL) 

CALL  SPURCYL(K»LL  »  IS ) 

NUMBFRsNUV,BER-l 

IF(NU^BER.Nf.C)  GO  TO  1211 

STOP 
END 
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SUBROUTINE  SPURCYL 

The  function  of  this  subprogram  is  to  perform  the  finite  difference 
calculations  for  the  set  of  diffusion-kinetics  equations  and  to  call  the 
subroutines  that  perform  the  interpolation  and  calculation  of  the  reaction 
integrals.   All  the  input  data  needed  in  this  program  are  brought  in  by 
the  use  of  a  COMMON  area  from  Program  DIFFUSE. 


Symbol 


Important  Program  Variables 


Explanation 


MARG(I,L) 


K 

NRADL, 
NRADLL 


TIME 

COEFFl(I), 
C0EFF2(I), 
C0EFF4(J,I), 
TIMEB(I) 


CZERO(I) 
CLEFT (I) 


Control  variable  determining  the  type  of  interpolation 
needed,  if  any 

Index  of  the  time  loop 


Indices  controlling  the  storage  location  of  the 
concentrations 

Current  reaction  time 


Factors  used  in  the  numerical  solution  of  the  finite 
difference  equations 

Initial  amount  of  radicals  of  species  I  present 

Amount  of  radicals  of  species  I  present  at  time  t 


The  subprograms  used  by  this  program  are: 

SUBROUTINE  INTMULT(I,L,M,LL) 
SUBROUTINE  INTGT1D(I ,LL,K,IO) 
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Logic  Diagram  for  Subroutine  SPURCYL 


c 


Enter  SUBROUTINE 
SPURCYL 


Determine 
MARG(I,L) 


K+l 


► — */yes\ 


CALCULATE 
TIME 


Calculate 
COEFFl(I), 
C0EFF2(I) 
C0EFF4(J,I) 
TIMEB(I) 


YES 


Calculate  new 
C(0,I) 


rCALL    INTMULI 


Calculate 
New  C(J,I) 


WRITE 
TIME, 
C(J,I) 


CALL  INTGT1D 


100 


101 


SUBRCU 

DIMENS 

DIMENS 

DIMENS 

DIMENS 

DIMENS 

COMMON 

COMMON 

COMMON 

COMMON 

COMMON 

l.S.NUM 

1Q.RADI 

1C0EFF4 

2  FORMAT 

3  FORMAT 
7  FORMAT 

DO  111 
N=NRXN 
DO  111 
IF( I.G 
I  I = IC  ( 
JJ=JC( 
IF( ITY 
IF(  ITY 
IF( ITY 
IF( ITY 
IF(D( I 
IF(D( I 
IFIDU 
GO  TO 


'L(K.LL.IS) 

)»D(5)»C(5CilO    )  *NRXNTERM( 5).IC(5.5).JC(5.5) 
) .A( 5) »  CZER0(5) .CLEFT ( 5) »RADIAL( 50) 

»WTR(5»50)»AXIAL(50)  »WN2( 50)  »MARG( 5*5 )  *CZLEFT(5  ) 


TINE  SPURCYI 
ION  RZER0(5 
ION  S ( 5 . 5  »  5 
ION  WNR(50 
ION  C0EFF1(5 
ION  SUM(5)  »' 
CZLOSS(5»5 


CMULT(5»50) »CZTLOSS( 5*5) »CLCSS(5>5) »DLTM(5 ) 
) 


39 


40 

1111 
173 


3432 
^433 


TAU(5 ) 

TIMEB( 5) 

TTIME 

SPUR»NUMRA 
ZPTS.SPURSE 
AL»WNR»WTR» 
»C0EFF5»TIM 
(44H0THE  TI 
(2X.8E16.8 ) 
(30H0THE  C 
1  I=1»NRADS 
TERM( I ) 
1  L=1.N 
T.NUMRAD.AN 
ItLJ 
IfL) 

PESOL.EQ.C. 
PESOL.EO.O. 
PESOL.ECO. 
PESOL.EO.O. 
)  .  N  E  .  D  (  I  I  )  . 
).EQ.D( JJ) . 
)  .NE.Dl JJ) . 
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D»NUMSCL.NUMRPTS»R2ERC»D»NRADS0L»C»NRXNTERM, IC.JC 
P»ALPHA.TINIT»DELTAZ.A»DELT4R,DELTAT»C2ERC.CLEFT» 
AXIAL.WNZ »PI  .  ITY PESOL» MARG »C0EFF1,CCFFF2. CCEFF3, 
E»DELTIMEtSUM»CMULTtC2TLCSS . 'LOSS ,DLTM ,C2L EFT 
ME-INCREASING  RADIAL  COEFFICIENTS  ARE) 

ONCENTRATIONS  OF  SPEC  I ES ♦ I  2 *4H  ARE*/) 


D. ITYPESOL.FQ.O)  GO  TO  173 


AND.  I  I .GT.NUMRAD) 
AND.  I  I .GT.NUMRAD) 
AND.JJ.GT.NL 4RAD) 
AND. JJ. GT.NUMRAD) 


Df II 1*01 1 1 

TAU(  I  I  )=TAU(  I  ) 
D( JJ)=D( I ) 
TAU( JJ)=TAU( I  ) 


CR.TAUU  )  .NE.TAUt  II  )  )  GO  TO  40 
AND.TAUI  I  ) .EQ.TAUI JJ)  )  MARG (  I  .L  )  =  1 
OR.TAU( I ) .NE.TAU( JJ) )  MARG ( I *L)=2 


IF(D(  I  ) .EO.D( JJ) .AND.TAU( I  )  .EQ.TAUI JJ) )  MARG(  I  ,L)=3 

IF(D(I).NE.D(JJ).OR.TAU(I).NE.TAU(JJ) )  MARG( I »L)=4 

CONTINUE 

LL  =  1 

DO  105  K=l  ,1300 

AN  =  K 

NA=(AN/10.-.0l ) 

NB=(AN/10.+.01 ) 

IF(NA.EO.NB)  GO  TO  3433 

FCRMAK20HCTHE  PROBLEM  TIME  IS»El0.4»SH  SECONDS.///) 

PRINT  6.  TIME 

DO  3432  I=1.NRADS0L 

PRINT  7. I 

PRINT  3.  (C( J. I  )  .  J=l  .NUMRPTS) 

CONTINUE 

CONTINUE 

NRADL=NRADSOL*LL 

NRADLL=NRADSOL*( 1-LL ) 

2  =  < 

TIME=TAU( 1)*(EXPF(2*DELTAT)-1.) 
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DC    IOC    I=1»NRADSCL 

IF( I.GT.NUMRAD.AND. ITYPESCL.EQ.O)    GC    TO    174 

TIMEB( I )=TIME+TAU( 1 ) 

CCEFF1 ( I )=DELTAT*TIMEB( I )/( ( DELTAR**2 ) * ( T IME+TAU ( I ) ) ) 

CCEFF2J  I  )=CCEFFK  I  )  *2  .*  (  1 .  +  ALPHA  ) 

DC  100  J=2tNUMRPTS 

2=J-1 

100  CCEFF4( J, I )=(0.5*ALPHA/Z+0.2  5*Z*DELTAR**2)*CCEFFK I ) 

174  CCNTINUE 
143  Z=K-1 

I  C=0 

DC  102  I=1»NRADSCL 

IF( I.GT.NUMRAD.AND. ITYPESCL.EQ.O)  GC  TC  175 

SUMM=0. 

N=NRXNTERM( I ) 

DC  101  J=l ,N 

ii=ic( i.j) 

JJ=JC( I .J) 

MI=II+NRADLL 
MJ=JJ+NRADLL 

101  SUMM   =SUMM   +S(  I  ♦  I  I  »JJ)*C(  1  »MI   )*Cd»MJ   ) 
MJ=I+NRADL 
MI=I+NRADLL 

102  C(  l.MJ)=C(  l.MI )+CCEFF2« I )* < C ( 2 »MI ) -C( 1 »MI ) )+T IMEB ( I ) *SUMM*DEL TAT 

175  CCNTINUE 

ic=o 

NR=NUMRPTS-1 

DC  104  I=1.NRADSCL 

IF( I.GT.NUMRAD.AND. ITYPESCL.EQ.O)  GC  TC  176 

N=NRXNTERM( I ) 

DC  106  L=1.N 

M=MARG( I .L) 
106  CALL  INTMULT( I .L»M,LL.l ) 

IF( I.GT.NUMRAD)  GC  TC  927 

CALL  INTGT1D( I .LL.K.IC) 
927  CCNTINUE 

MI=I+NRADLL 

MJ=I+NRADL 

DC  104  J=2.NR 

SUMM=0. 

DC  115  L=1.N 
115  SUMM   =SUMM   +CMULT(L.J) 

104  C(J»MJ)=C(J.MI  )+CCEFFl (  I  ) * ( C ( J+l »M I ) ~2 . *C ( J »M I  ) +C ( J~l »M I )  ) 
1+CCEFF4( J, I )*(C( J+l »MI )-C( J-l »MI ) )+TlMEB( I ) *SUMM*DELTAT 

176  CCNTINUE 
IF(IC.EQ.I)  RETURN 
Z  =  0. 

DC  120  I=1,NUMRAD 

IF(CLEFT( I )/CZERC( D.LT.0.02)  Z^Z+1. 
120  CCNTINUE 
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Y=NUMRAD 

IF(Z.EQ.Y)  RETURN 
IF(Ll.EQ.l)  GC  TO  330 

LL=1 

GC  TO  105 
330  LL=0 

IF(TIME.GE.TTIME)  RETURN 

105  CONTINUE 
RETURN 
END 
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SUBROUTINE  INTMULT(I,L,M,LL) 

The  function  of  this  subprogram  is  to  perform  the  interpolation  and 
multiplication  necessary  for  calculation  of  the  reaction  terms  in  the 
finite  difference  equations  and  for  the  calculation  of  the  reaction 
integrals. 

The  interpolation  scheme  used  is  the  second  order  forward  Gregory- 
Newton  interpolation  formula  given  as 

f  (x)  =   f  (x0  +  rh)  =  f0  +  rAf0  +  -^-^-A2f0  (80) 

2! 

where:    x  =  the  tabulated  value  of  position  just  below  the  interpolation 

point 

h  =  the  magnitude  of  the  grid  spacing 

r  =  a  number  between  zero  and  one  such  that  x  +  rh  is  the 

o 

value  of  the  interpolation  point 
AfQ  =  first  forward  difference 
A2fQ  =  second  forward  difference 

Important  Program  Variables 


Symbol 


Explanation 


M  Control  variable  equal  to  MARG(I,L)  which  determines 

the  type  of  interpolation  needed 

H  Number  corresponding  to  r  in  Eq.  (80) 

CA,  CB         Values  of  the  interpolated  concentrations 

CMULT(I,L)      Value  of  the  interpolated  reaction  term 


Logic  Diagram  for  Subroutine  INTMULT 
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c 


Enter 
Sub.  Intmult 


Interpolate  to 
determlne  CA 


Set 
CA=C(j,N) 


Set 
CB=C(J,n) 


Interpolate  to 
determine  CB 


Calculate 
CMULT(L.J) 


J  =  J  +  1 


/  NO  V < 
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SUBROUTINE  INTMULK I »L»M»LL»JZ) 

DIMENSION  R2ERO(5)  »0 ( 5 ) »C( 50 » 10    )  »NRXNTERM( 5 ) » IC ( 5 ♦ 5 ) ♦ JC ( 5 , 5 ) 

DIMENSION  S( 5»5»5)  »A(5)  »  CZERC ( 5 ) »CLEFT ( 5 )  .RADI AL ( 50 ) 

DIMENSION  WNR(50) *WTR( 5.50) » AXIAL ( 50) »WNZ( 50 ) »MARG ( 5 » 5 ) » CZLEFT ( 5 ) 

DIMENSION  C0EFFK5)  »COEFF2(5)  >C0EFF3(5)  »C0EFF4  (  50  »  5  )  »C0EFF5  (5*50) 

DIMENSION  SUM (5 ) ,CMULT(5»50) >CZTLOSS( 5»5) ,CLOSS(5>5) »DLTM<5  > 

COMMON  CZLOSS(5,5) 

COMMON  TAU(5) 

"COMMON  TIMEB(5)  '  '  "   ': •••">:' 

COMMON  TTIME 

COMMON    SPURtNUMRAD»NUMSCL»NUMPTS»RZERC»D»NRADSCL»C»NRXNTERM» IC»JC 
ltS»NL'MZPTS.SPURSEP.ALPHA»TINIT»DELTAZ»A»DELTAR»DELTATtCZERCfCLEFTt 
lQtRADIAL»WNRtWTRtAXIAL»WNZ»Pl  *  I TYPESOl ,MARG »C0EFF1 , C0EFF2 »C0EFF3 * 
1CCEFF4»CCEFF5»TIME»DELTIME»SUM»CMULT»CZTLCSS,CLCSS»DLTM,CZLEFT 
2    FORMAK  15  ) 

NR=NUMRPTS-1 

I  I  =  IC( I »L) 

JJ=JC( I »M 

LI=NRADSOL*( 1-LL) 

M I  =  I  I +L  I 

MJ=JJ+LI 

R=(  (D( I >* (TAU(  I )+TIME) / (D(  II  )*(TAU( I  I  )+TlME) ) ) **0. 5 ) *DELTAR 

RR=( (DU )*(TAU( I )+TlME)/(D( JJ)*<TAU( JJJ+TIME) ) ) **0. 5 ) *DELTAR 

DO  ICO  J=2»NR 

CA  =  0. 

CB  =  0. 
GO  TO  (302*301, 300»300) »M 

300  Z=J-1 
RINTERP=Z*R 

N=(RINTERP/DELTAR     ) 
IF(N.GT.NUMRPTS)  GO  TO  307 
Z  =  N 

Z=Z*DELTAR 

H=(RINTERP-Z)/DELTAR 

N  =  N  +  1 

CA=H*(2.-H)*C(N+1»MI    ) +0 • 5* ( H-l . )* ( H*C ( N+2 »M I    )+(H-2.)* 
1C(N»MI    )) 

GO  TO  305 
307  CA=C(NUMRPTS»MI ) 
305  IF(M.E0.3)  GO  TO  303 

301  Z  =  J-1 
RINTERP=Z*RR 
N=(RINTERP/DELTAR     ) 
IF(N.GT.NUMRPTS)  GO  TO  308 
Z  =  N 

Z=Z*DELTAR 

H=(RINTERP-Z)/DELTAR 
N  =  N  +  1 

CB  =  H*(2.-H)*C(N  +  1»MJ    )+0 . 5* ( H-l .  )  * ( H*C ( N  +  2 *MJ    )+(H-2.)* 
1C(N,MJ) ) 
GO  TO  306 
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308  CB  =  C(NUMRPTS»V!J) 

306  IF(M.E0.2)  GO  TC  304 

CMULT(L*J)=5( [ill »JJ)*CA*CB 

GO  TO  100 

303  CMULT(L*J)=S( I » II »JJ)*CA*C< JtMJ    ) 
GO  TC  100 

304  CMULT(L»J)=S( I»II»JJ)*C(J»MI    )*CB 
GO  TC  100 

302  CMULT(L»J )=S(  I »I  I » JJ)*C( J»MI    )*C(J»MJ 
100  CONTINUE 
RETURN 
END 
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SUBROUTINE  INTGT1D(I ,LL,K,IO) 

The  function  of  this  subroutine  is  to  calculate  the  reaction  loss 
integrals  and  the  current  amounts  of  free  radicals  present.   These 
results  are  then  printed  out  at  every  tenth  time  step.  The  volume 
integrals  are  determined  numerically  by  a  Simpson's  Integration  scheme 
with  weights  given  by  the  FATES  subprogram. 

Important  Program  Variables 


Symbol 


Explanation 


CL0SS(I,L)      Cumulative  amount  of  the  I   species  created  by  the 
Ltn  reaction  term 

CTOTAL         The  first  value  determined  is  that  for  the  amount  of 
the  I*-"  species  present  at  time  t.   The  second  value 
determined  is  that  for  the  sum  of  the  amount  of  the 
I   species  present  and  the  total  amount  lost  by 
reaction.   This  value  should  remain  close  to  the 
value  for  the  initial  amount  of  the  1^   species  pre- 
sent.  If  too  large  a  deviation  develops  this  sub- 
program will  cause  the  main  program  to  terminate  this 
particular  calculation. 


Logic  Diagram  for  Subroutine  INTGT1D 
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c 


ENTER 
SUB.  INTGT1D 


Calculate 
CL0SS(I,L) 


WRITE 
CLOSS(I.L) 
CTOTAL 


Calculate 
CTOTAL 


SET        NRXNTERM(I) 

CTOTAL=CTOTAL  -[CL0SS(I,L) 
L=l 


WRITE 
CTOTAL 


Return 
.SUB,  INTGT11 


SET 
10  =  1 
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DIMENSION 

DIMENSION 

DIMENSION 

DIMENSION 

COMMON 

COMMON 

COMMON 

COMMON 

COMMON 


SUBROUTINE  I NTGT ID ( I »LL »K * 
DIMENSION  R2ER0(5)  »D<5)  »C( 
S(5»5»5) »A(5)  » 
WNR( 50) »WTR( 5*50 
CCEFFK5)  >C0EFF2 
SUM(5) »CMULT(5»5 
CZLOSS(5,5) 
TAU(5  ) 
TIMEB(5) 
TTIME 

SPUR,NUMRAD,NUMSOL, 
ltS»NUMZPTStSPURSEP» ALPHA »T 
1Q »R AD  I AL , WNR ♦ WTR  » A  X  I AL , WNZ 
lCCEFF4»CCEFF5fTIME»DELTIME 

3  FORMAT(2X,8E16.8) 

4  F0RMAT(24HCTHE  REACTION  LO 

5  FORMAT (35H0THE  CURRENT  NUM 

6  FORMATdlH  THE  SUM  IS»E15. 
I I=I+NRADSOL*( 1-LL) 

AN  =  K 

NA=(AN/10.+.01) 

NB=(AN/10.-.0l) 

Z=4.*PI*TIMEB( I )*( (D( I )*(T 

N=NRXNTERM( I ) 

DO  102  L=1»N 

SUMM=0. 

DO  100  J=2»NUVRPTS 

Y  =  J-1 

10^  SUMM=SUMM+CMULT(L» J)*WNR( J 

102  CLOSS( I »L)=CLOSS( I »L)+SUMM 
. IF(NA.EQ.NB)  GO  TO  301 

Z=Z/(TIMEB( I )*DELTAT) 

SUMM=0. 

DO  101  J=2»NUMRPTS 

Y  =  J-1 

101  SUMM=SUMM+C( Jt II )*WNR( J)*( 
CTCTAL=Z*SUMM 

z=o. 

DO  103  L=1»N 

103  Z=Z+CL05S( I iL) 
PRINT  4 

PRINT  3»(CL0SS( I »L)  ,L=1,N) 
302  PRINT  5»CT0TAL 

CTOTAL=CTOTAL-Z 

PRINT  6>CT0TAL 

IF(ABSF(CTOTAL/CZERO( I )-l. 
301  RETURN 
END 


10) 

50,  10 


)  * 

cz 

)  *  AXIAL ( 5C 


NRXNTERI 
ER0(5) » 
)  ,WNZ(  5' 


M(5) t IC( 5»5) * JC(5»5) 
CLEFT(5) »RADIAL( 50) 
0)  »MARG( 5»5 )  »CZLEFT( 5) 
(5)*C0EFF3(5) ,C0EFF4 ( 50 » 5 ) »C0EFF5 (5*50) 
0) »CZTLOSS( 5,5  J »CL0SS(5»5) *DLTM(5 ) 


NUMRPTStRZ 
INIT, DELTA 
,PI ,ITYPES 
,SUM,CMULT 


ERO,D,NRADSOL,C,NRXNTERM,  ICJC 


SSES 

BER  : 
8) 


Z,A,DEL 
OL,MARG 
,CZTLOS 

) 


T AR, DE L T AT, CZERO, CLEFT, 

»C0EFF1,C0EFF2,C0EFF3, 

S,CLOSS,DLTM,CZLEFT 


ARE,/ 
IF  PARTICLES  IS,E15.8,/) 


AU( I  J+TIME)  )**1.5)*DELTAT*DELTAR**2 


)*(Y**2) 
*Z 


Y**2) 


).GT.0.30)  10=1 
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SUBROUTINE  WATES (L,NUMRPTS , RADIAL, WNR) 

The  function  of  this  subprogram  is  to  generate  the  necessary  weights 
for  a  Simpson-type  numerical  integration  for  either  a  linear  or  a  logar- 
ithmic scale.   If  the  number  of  points  is  even,  a  cubic  approximation  is 
used  to  subtract  out  the  area  of  the  center  portion  of  the  curve  since 
the  Simpson  integration  counts  this  area  twice.   If  there  are  only  two 
points,  the  area  is  determined  using  the  trapezoidal  rule.   If  there  is 
only  one  point  the  area  is  set  to  zero. 

Important  Program  Variables 


Symbol 


Explanation 


NUMRPTS 
RADIAL (J) 

WNR (J) 


Control  integer 

■  1  for  a  linear  scale 

>  1  for  a  logarithmic  scale 

Number  of  points  considered  for  the  integration 

Variable  which  locates  the  abscissa  points  for  the 
integration 

Value  of  the  weights  determined  by  this  subprogram 
to  be  used  in  the  numerical  integration 


Logic  Diagram  for  Subroutine  WATES 
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Enter  Sub. 
WATES 


f\ 


RETURN 
ISub.  Wates 


Calculate 
factOi.'  between 
Points 


Multiply 
WNR(J)  by 
Proper  factor 


Calculate 
Interval  between 
Points 


SET 
Interval= I  Interval 


Calculate 
WNR(J) 


*— -f  NO  J 


Calculate 
WNR(l)  and  WNR(2) 


Multiply 
:NR(1)  and  17NR(2) 

bv  proper  factor 


113 


SUBROUTINE  FATESt IWT »NWT »WTAB »WATES ) 
DIMENSION  WTAB(50) »WATES(50) 
819  WTA=NWT 

IF(NWT-2.GE.C)  GO  TO  39 
19  WATESd  )  =  C. 

GO  TO  259 
39  IF( IWT-2.GE.0)  GO  TO  79 
59  WTDEL=(WTAB(  1  )-WTAB(NWT)  )/(WTA-l. ) 

GO  TO  99 
79  IF( IWT.EQ.3)  GO  TO  3333 

WTDEL=ALCG(WTAB( 1 )/WTAB(NWT) )/(WTA~l. ) 
99  IF(WTDEL.GE.O. )  GO  TO  990 
119  WTDEL=~WTDEL 
GO  TO  990 
^33  3  WTDEL=EXP(WTAB( 1)/WTAB(NWT) )/(WTA~l. ) 

990  IF(NWT-2)  259,1190,139 
1190  WATESl 1 )=  .5*WTDEL 
WATES(2)=WATES(1) 
GO  TO  199 
139  NWTA=(WTA/2.+»ll 
NWTB=(WTA/2.-.l ) 
NWTC=(WTA/4.+.l ) 
NWTD=(WTAM.-.l  )  ■ 
WATES( 1 )=WTDEL/3. 
WTC  =  WATES(  1) 
WATES(NWT)=WATES( 1  I 
DO  159  I=1»NWT3 
WATES( 1+1 )=WTDEL+WTC 
INDX=NWT-I 

WATES( INDX)=WTDEL+WTC 
159  WTC=-WTC 

WTD=l./24. 

IF(NWTC-NWTD.LE.O)  GO  TO  1790 
1590  WTO=-WTD 

1790  IF(NWTA-NWTB.LE.O)  GO  TO  199 
179  WATES(NWT3)=WATES(NWTB)-WTD*WTDEL 

WATES(NWT3+1 ) =WATES( NWTB+1 ) + 5 . *WTD*WTDEL 
WATES(NWTB+3  )  =WATES<  NWTB) 
WAT ES(NWTB  +  2)=a,ATES(  NWTB  +  1  ) 
199  IFCIWT-2.LT. Oj  GO  TO  259 
21°  DO  239  1  =  1  ,NWT 
239  WATESl I )=WATE5( I 1*WTAB( I ) 
259  RETURN 
END 
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A. 2.1  Sample  Input  Data  for  Program  Diffuse 


12 

.10000000E-04 

1     3    30 

1 
.70700000E-07 

.2OOOOOOOP-C4   .20000000E-04   • ?C0C0000e-04   .2OOOOOOOE-O4 
3     2     2     1 

1  4 
..10000000E-10 

2  3 
.10000000E-10 

1     1 

..lOOOOCOOE-10 

1  4 

.10000000E-10 

2  3 
-.lOOOOOOOE-10 

1  1 
•5C000000E-11 

2  3 
-.lOOOOOOOE-10 

1  4 

■•lOOOOOOOE-10 
.60000000E+01       .O00000O0E+00       . 60235000E+19       . 60235C00E+19 


Note:      Each  line   represents   one  computer   card. 
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A. 3   Selected  Data  for  Figures  1-12 


Run 

1 

Time    (sec) 

T 
NRS(t) 

T 

2NR2(t) 

yo 

h 

1.561xlO"U 

.006 



.345 

5.771 

6.122 

3. 525x10" U 

.012 



.620 

5.456 

6.088 

5.764xl0~U 

.019 



.846 

5.193 

6.058 

8.316xlO~H 

.026 



1.035 

4.970 

6.031 

1.122xl0~10 

.034 



1.197 

4.776 

6.006 

2.262xl0~10 

.062 

.002 

1.566 

4.314 

5.940 

A.676xl0"10 

.115 

.005 

1.887 

3.867 

5.865 

6. 448x10" 10 

.152 

.007 

2.002 

3.684 

5.831 

1.015xl0"9 

.226 

.013 

2.134 

3.437 

5.784 

2.067xlO"9 

.420 

.031 

2.280 

3.048 

5.712 

4.090xl0~9 

.755 

.106 

2.364 

2.635 

5.648 

6.988xl0"9 

1.174 

.250 

2.402 

2.275 

5.601 

l.OAOxlO"8 

1.610 

.467 

2.420 

2.007 

5.570 

2.013xl0~8 

2.663 

1.232 

2.436 

1.655 

5.522 

4.427xl0"8 

4.922 

3.376 

2.446 

1.485 

5.477 

6.560xl0"8 

6.834 

5.284 

2.449 

1.458 

5.458 

1.108xl0-7 

10.818 

9.268 

2.453 

1.431 

5.433 

2.131xl0"7 

19.672 

18.123 

2.456 

1.399 

5.404 

A.lOOxlO"7 

36.345 

34.798 

2.458 

1.369 

5.374 

6.920xl0~7 

59.784 

58.237 

2.459 

1.345 

5.352 
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Run   2 


Time    (sec) 

NRS(t) 

T 
NpCt) 

2NR    (t) 
R2 

yo 

I* 

1.561xlO~U 

.061 



.343 

5.723 

6.126 

3. 525x10" U 

.123 

.002 

.610 

5.359 

6.091 

5.764xl0"U 

.188 

.004 

.826 

5.051 

6.061 

8. 316x10" H 

.257 

.009 

1.004 

4.781 

6.033 

1.122xl0"10 

.330 

.015 

1.153 

4.540 

6.008 

2.262xl0"10 

.590 

.052 

1.478 

3.928 

5.943 

4.676xl0~10 

1.060 

.170 

1.731 

3.252 

5.872 

6.448xl0"10 

1.936 

.569 

1.892 

2.539 

5.798 

1.015xl0"9 

3.301 

1.574 

1.966 

2.044 

5.737 

2.067xl0~9 

5.583 

3.730 

2.006 

1.821 

5.689 

4.090xl0~9 

8.699 

6.838 

2.026 

1.771 

5.658 

6.988xl0~9 

12.313 

10.456 

2.038 

1.740 

5.636 

1.040xl0"8 

22.412 

20.561 

2.053 

1.695 

5.599 

2.013xl0"8 

46.907 

45.060 

2.065 

1.646 

5.557 

4.427xl0-8 

68.173 

66.329 

2.069 

1.624 

5.536 
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Run   3 


Time    (sec) 

<s(t) 

Np(t) 

2NL(t) 

yo 

z. 

1.561X10"11 

.586 

.029 

.318 

5.284 

6.159 

3.525xl0~n 

1.137 

.126 

.529 

4.567 

6.107 

5.764xlO-11 

1.673 

.294 

.677 

4.010 

6.067 

8.316xlO~U 

2.206 

.536 

.786 

3.576 

6.033 

1.122xlO~10 

2.751 

.855 

.869 

3.240 

6.005 

2. 262x10" 10 

4.565 

2.306 

1.032 

2.650 

5.942 

A.676xl0~10 

7.924 

5.579 

1.163 

2.380 

5.889 

6.448xl0~10 

10.030 

7.960 

1.211 

2.320 

5.868 

1.015xl0~9 

15.165 

12.848 

1.270 

2.252 

5.838 

2.067xl0~9 

28.725 

26.438 

1.341 

2.164 

5.791 

4.090xl0~9 

54.227 

51.963 

1.390 

2.093 

5.746 

6.988xl0"9 

90.068 

87.817 

1.417 

2.043 

5.711 
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Run  4 


Time  (sec) 

NRS(t) 

T 
Np(t) 

2NR2(t) 

Vt) 

k 

1. 561x10" H 

4.440 

1.810 

.193 

3.304 

6.127 

3. 525x10" U 

7.938 

5.116 

.278 

2.938 

6.038 

5.764xl0"H 

11.678 

8.874 

.350 

2.866 

6.020 

8. 316x10" H 

15.819 

13.044 

.414 

2.812 

6.007 

1. 122x10" 10 

20.434 

17.686 

.473 

2.774 

5.994 

2. 262x10" 10 

37.899 

35.221 

.620 

2.661 

5.959 

4. 676x10" 10 

73.428 

70.819 

.763 

2.841 

5.914 

6. 448x10" 10 

98.926 

96.343 

.819 

2.491 

5.892 

1.015xl0"9 

151.336 

148.787 

.887 

2.425 

5.861 
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Run   5 


Time  (sec) 

NRS(t) 

T 
Np(t) 

2NR2(t) 

yo 

h 

1.561x10" 

.586 

.003 

.317 

5 

.264 

6.163 

3.525xl0"U 

1.130 

.016 

.524 

4 

.477 

6.115 

5.764xl0~U 

1.647 

.039 

.662 

3 

.806 

6.077 

8.316xlO"U 

2.141 

.074 

.756 

3 

.221 

6.044 

-10 
1.122x10 

2.614 

.122 

.819 

2 

.706 

6.017 

2.262xl0~10 

3.891 

.355 

.906 

1 

.510 

5.952 

-10 
4.676x10 

5.240 

.939 

.931 

.668 

5.900 

-10 
6.448x10 

5.831 

1.385 

.934 

.506 

5.886 

1.015x10" 

6.828 

2.317 

.936 

.427 

5.875 

2.067x10" 

9.453 

4.935 

.939 

.408 

5.865 

4.090x10" 

14.432 

9.918 

.940 

.402 

5.857 

6.988x10" 

21.511 

17.000 

.941 

.398 

5.850 

1.040x10" 

29.806 

25.298 

.942 

.396 

5.845 
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Run  6 


Time    (sec) 

NRS(t) 

T 
Np(t) 

2NR2(t> 

NR(t) 

^R 

1.561xlO~H 

.061 

.003 

.343 

5.725 

6.126 

3.525xl0"n 

.123 

.013 

.611 

5.369 

6.090 

5.764xl0~n 

.188 

.032 

.828 

5.075 

6.059 

8.316xl0-11 

.258 

.061 

1.009 

4.826 

6.031 

1.122xl0~10 

.332 

.100 

1.161 

4.613 

6.006 

2.262xlO~10 

.600 

.296 

1.505 

4.131 

5.940 

4. 676x10" 10 

1.114 

.779 

1.804 

3.727 

5.866 

6.448xl0"10 

1.473 

1.139 

1.913 

3.586 

5.834 

1.015xl0"9 

2.201 

1.872 

2.042 

3.417 

5.788 

2.067xl0"9 

4.189 

3.868 

2.195 

3.201 

5.717 

4.090xl0"9 

7.869 

7.553 

2.295 

3.039 

5.650 

6.988xl0~9 

12.988 

12.676 

2.351 

2.936 

5.600 

1.040xl0~8 

18.905 

18.594 

2.383 

2.869 

5.562 
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Run  7 


Time    (sec) 

NRS(t) 

T 
Np(t) 

2NR    (t) 
K2 

yo 

h 

1. 561x10" n 

.089 



.743 

8.372 

9.203 

3. 525x1 0"11 

.177 

.003 

1.285 

7.684 

9.143 

5.764xl0~n 

.266 

.008 

1.704 

7.132 

9.093 

8.316xl0"u 

.358 

.015 

2.038 

6.670 

9.051 

1.122xlO-10 

.456 

.026 

2.310 

6.273 

9.014 

2.262xlO~10 

.796 

.082 

2.887 

5.319 

8.921 

4.676xlO-10 

1.406 

.251 

3.327 

4.341 

8.823 

6. 448x10" 10 

1.802 

.406 

3.463 

3.921 

8.781 

1.015xl0~9 

2.541 

.785 

3.603 

3.364 

8.724 

2.067xlO"9 

4.314 

2 

.092 

3.731 

2.690 

8.644 

4.090xl0~9 

7.271 

4 

.878 

3.798 

2.389 

8.581 

6.988xl0~9 

11.296 

8 

.890 

3.833 

2.300 

8.539 

1.040xl0~8 

15.959 

13 

.558 

3.853 

2.256 

8.510 

2.013xl0~8 

28.987 

26 

.594 

3.877 

2.194 

8.463 

4.427xlO~8 

60.579 

58 

.195 

3.897 

2.128 

8.409 
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Run  8 


Time    (sec) 

NRS(£) 

T 
Np(t) 

2NR2(t) 

yo 

k 

1. 561x10" U 

.116 



1.274 

10.897 

12.286 

3. 525x10" n 

.226 

.005 

2.149 

9.825 

12.196 

5.764xl0"n 

.335 

.012 

2.799 

9.003 

12.126 

8. 316x10" n 

.448 

.023 

3.305 

8.339 

12.068 

1.122xl0"10 

.565 

.039 

3.711 

7.782 

12.019 

2. 262x10" 10 

.970 

.117 

4.550 

6.497 

11.900 

4.676xl0"10 

1.691 

.334 

5.176 

5.246 

11.779 

6. 448x10" 10 

2.158 

.525 

5.370 

4.726 

11.728 

1.015xl0~9 

3.032 

.983 

5.568 

4.043 

11.660 

2.067xl0"9 

5.128 

2.532 

5.751 

3.216 

11.564 

4.090xl0~9 

9.592 

6.772 

5.858 

2.796 

11.475 

6.988xl0"9 

13.353 

10.524 

5.893 

2.717 

11.439 

1.040xl0"8 

21.171 

18.350 

5.928 

2.644 

11.393 

2.013xl0~8 

34.133 

31.321 

5.954 

2.583 

11.349 

4.427xl0~8 

71.231 

68.430 

5.981 

2.503 

11.285 

6.560xl0"8 

103.438 

100.642 

5.991 

2.466 

11.254 
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Run  9 


Time  (sec) 

slsM 

T 
Np(t) 

2NR  (C) 

yo 

[R 

1.561x10" 

.063 



2 
.004 

6.036 

6.102 

3.525xl0~U 

.129 

.001 

.007 

5.941 

6.076 

5.764xlO-11 

.202 

.004 

.010 

5.841 

6.050 

8.316xlO"U 

.282 

.007 

.013 

5.738 

6.025 

1.122xlO~10 

.369 

.013 

.016 

5.629 

6.001 

2.262xlO"10 

.693 

.048 

.022 

5.262 

5.929 

4. 676x10" 10 

1.323 

.180 

.027 

4.668 

5.838 

6. 44 8x10" 10 

1.750 

.317 

.029 

4.333 

5.795 

1.015xl0"9 

2.570 

.785 

.031 

3.817 

5.734 

2.067xl0~9 

4.585 

2.078 

.033 

3.102 

5.643 

4.090xl0~9 

7.982 

5.218 

.034 

2.771 

5.569 

6.988xl0~9 

12.644 

9.843 

.035 

2.685 

5.520 

1.040xl0~8 

18.079 

15.272 

.035 

2.644 

5.486 

2.013xl0"8 

33.354 

30.543 

.035 

2.585 

5.431 

4.427xl0~8 

70.613 

67.801 

.036 

2.519 

5.367 

6.560xl0~8 

103.066 

100.254 

.036 

2.487 

5.335 
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Run  10 


Time    (sec) 

NRS(t) 

T 
Np(t) 

T 

2nrF 

yt) 

Ir 

1.561xl0-11 





.035 

5.776 

6.122 

3.525xl0"n 

.001 



.621 

5.465 

6.088 

5.764xl0"n 

.002 



.848 

5.208 

6.058 

8.316xl0~n 

.003 



1.038 

4.990 

6.031 

1. 122x10" 10 

.004 



1.201 

4.801 

6.006 

2.262xl0~10 

.007 



1.575 

4.358 

5.939 

A. 676x10" 10 

.013 

.002 

1.906 

3.948 

5.864 

6. 448x10" 10 

.017 

.004 

2.025 

3.791 

5.829 

1.015xl0"9 

.026 

.007 

2.167 

3.596 

5.781 

2.067xl0"9 

.048 

.022 

2.332 

3.348 

5.706 

4.090xl0~9 

.088 

.057 

2.440 

3.166 

5.637 

6.988xl0"9 

.143 

.111 

2.500 

3.052 

5.584 

1.040xl0"8 

.206 

.175 

2.534 

2.980 

5.546 

2.013xl0~8 

.381 

.350 

2.576 

2.877 

5.483 

4.4x7xl0~8 

.798 

.767 

2.609 

2.773 

5.412 

6.560xl0"8 

1.158 

1.127 

2.602 

2.726 

5.377 

1.108xl0"7 

1.906 

1.875 

2.632 

2.669 

5.332 

2.131xl0~7 

3.561 

3.530 

2.643 

2.603 

5.277 
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Run  11 


Time    (sec) 

NRS(t) 

T 
Np(t) 

2NR2(t) 

yo 

k 

1.561xlO~H 

.061 

.025 

.343 

5.743 

6.122 

3.525xlO-11 

.124 

.083 

.616 

5.430 

6.087 

5.764xlO~n 

.190 

.150 

.840 

5.177 

6.057 

8.316xlO-11 

.261 

.222 

1.029 

4.962 

6.030 

1.122xlO~10 

.338 

.300 

1.191 

4.777 

6.005 

2.262xlO~10 

.617 

.580 

1.563 

4.430 

5.940 

4.676xlO"*10 

1.158 

1.121 

1.892 

3.936 

5.864 

6. 4 4 8x10" 10 

1.536 

1.500 

2.011 

3.783 

5.830 

1.015xl0~9 

2.300 

2.263 

2.153 

3.593 

5.782 

2.067xlO"9 

4.380 

4.345 

2.319 

3.354 

5.707 
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Run   12 


Time  (sec) 

NRS(t) 

T 
Np(t) 

2NR2(t) 

yo 

k 

8.377x10 

.284 



.020 

5.728 

6.032 

1.775x10 

.489 



.039 

5.497 

6.024 

2.719x10 

.647 



.056 

5.317 

6.020 

4.629x10 

.882 

.002 

.089 

5.044 

6.014 

7.548x10 

1.127 

.003 

.134 

4.750 

6.008 

1.053x10 

1.306 

.006 

.177 

4.527 

6.004 

2.097x10 

1.700 

.016 

.302 

4.009 

5.995 

4.060x10 

2.117 

.040 

.474 

3.433 

5.984 

6.997x10 

2.537 

.084 

.640 

2.878 

5.972 

1.046x10 

2.942 

.144 

.761 

2.402 

5.961 

2.028x10 

3.882 

.344 

.914 

1.490 

5.942 
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ABSTRACT 

The  objectives  of  this  work  were  two-fold.   The  first  was  to 
develop  a  general  computational  scheme  applicable  to  the  solution  of  the 
equations  representing  the  diffusional  and  kinetic  behavior  of  reacting 
species  in  a  single  cluster  along  the  track  of  a  particle  of  ionizing 
radiation.   The  second  was  to  apply  this  method  to  the  investigation  of 
the  spatial  and  time  behavior  of  the  concentrations  of  reactants  in  a 
typical  radiation- induced  chemical  chain  reaction.   This  investigation 
was  done  for  a  variety  of  initial  concentrations  and  chemical  reaction 
rate  constant  parameters. 

The  extent  of  reaction  as  a  function  of  time  for  each  reacting  species 
was  determined  for  each  set  of  parameters.   This  enabled  general  conclu- 
sions to  be  drawn  for  the  effects  of  the  parameters  on  the  possibility 
of  cluster  overlapping  in  a  chain  reaction.   These  data  were  also  used 
to  qualitatively  determine  the  effects  of  variation  of  dose  rate  and  to 
give  an  indication  of  the  dose  rates  at  which  spur  overlap  may  become 
important. 

The  spatial  and  time  dependencies  of  concentration  were  determined 
for  one  set  of  parameters  in  order  to  indicate  the  feasibility  of  utilizing 
a  modified  form  of  the  prescribed-diffusion  hypothesis. 


